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ABSTRACT 



Digital control offers exciting possibilities for ac- 
curacy and flexibility that cannot be achieved by analog 
means. The USQ-2Q General Purpose Digital Computer is cn 
board several ships today for the Navy Tactical Data System, 
and an effort is made to show how such a computer might also 
be employed directly in closed-loop shipboard fixe control 
systems. Some of the theory of sampled data systems and 
digital compensation is discussed. A digital computer pro- 
gram which simulates the action of a sampled data system 
with a digital controller is employed to experimentally de- 
termine the behaviour of various physical plants subject to 
digital control. An evaluation is made of the noise rejec- 
tion capabilities, steady state error, and transient control 
of a digital process. A full description and instructions 
for its use is given for the simulator program, which has 
proved to be very helpful in making theoretical studies of 
control system behaviour. 

The writer wishes to express his appreciation for the 
assistance and encouragement given him by Professor Mitchell 
L. Cotton of the U„ S„ Naval Postgraduate School and John B. 
Slaughter of the U„ S c Navy Electronics Laboratory in this 
investigation. 
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Introduction 



A control system which contains a digital computer 
may be classed as a sampled-data system* Such systems are 
becoming increasingly more common, both in military opera- 
tions and in the commercial field. By the use of sampling 
in control systems it is possible to build simple, sensi- 
tive, and efficient power control devices. Sampling per- 
mits the control of tremendous power by sensitive control 
elements without excessive amplification, and also mini- 
mizes the loading effect upon sensitive instruments. 
Utilizing sampled data and digital components in a control 
system allows time sharing of important parts of the system. 
This is a very significant advantage of digital control 
systems. ' /j ' / 

Since sampled-data can readily be coded, the data 
signals in digital control systems are received and trans- 
mitted in pulse-code form, and they will provide almost 
error-free channels for transmission through noisy media. 
Consequently, the only noise in the transmission of pulse- 
coded signals is the quantization error. Likewise, a digi- 
tal computer in a control system allows data processing of 
control information so that the flexibility and versatility 
of the digital computer can be utilized to improve the per- 
formance of the control system. Furthermore, digital con- 
trol techniques make feasible the system compensation by 
non-linear programming and adaptive, or self -optimizing, 
control. This is another significant advantage achievable 
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with digital control 



/ 2 / 



The principal features of a typical system containing 



a digital computer are shown in block diagram form in Fig- 
ure 1 . 
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Figure 1. Typical block diagram of a hybrid control system. 

A characteristic feature of this system is the fact that 
data will appear at one or more places within it as a se- 
quence of numbers or fixed signal levels that do not change 
appreciably except at fixed instants in time. The process 
that converts continuous data into such a sequence is called 
the sampling process, and the mechanism that implements the 
process is called an analog- to-digital converter (ADC). The 
process that converts the output number sequence of the digi- 
tal computer to a more or less continuous analog quantity 
is called the desampling process, and the mechanism accomp- 
lishing this process is called the digital-to-analog con- 
verter (dac ) 

The quantities T^, T 2 , and T^ are the periods of the 
samplers shown. In almost every case T^ is the same as T 2 , 
but the relation between T^ and T^ may be chosen by the de- 
signer. The simplest relation, of course, is when T ? equals 
T^. The system is then called single - rate since all discrete 
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number sequences appearing within the system have the same 



sampling period. A multi - rate system occurs when T ? and 
are not the same. Analysis is simplified when and 
are related by an integer and we have 11 synchronous" 
sampling. It is also simpler to analyze a multi-rate sys- 
tem wherein T 3 is less than T 2 (ie., T 3 - T^/n n = 2,3,...) 
which is fortunate, since most of the advantages of a multi- 
rate system are derived when the output of the digital con- 
troller occurs at a faster rate than the input. 

The experimental results described in this thesis are 
obtained from a digital computer program that simulates a 
hybrid control system from beginning to end. See Appendix I 
for the description of this program. In an actual control 
system the only operations that would perforce be external 
to the digital computer would be the plant itself and the 
ADC. In the investigation conducted by John Slaughter at 
the U. S. Navy Electronics Laboratory, San Diego, California, 
a Packard-Bell Model M2 analog- to-digital converter is em- 
ployed, and a 100 micro-second conversion and settling time 
is allowed. For a discussion of a proposed analog- to-digital 
system for use with the CDC 1604 digital computer see Ref- 
erence 6. 

The generation of the error signal by comparison be- 
tween input and output signals can be accomplished either 
externally to the digital computer or by the digital com- 
puter itself. The latter case would require a two channel 
ADC and a multiplexing system for sampling respectively the 
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input and the output signals. 

Data reconstruction by the DAC can also be performed 
by an external device and can be of various degrees of com- 
plexity. However, the fact that the plant in control systems 
is usually low-pass in frequency response makes the use of 
overly complex data-reconstruction systems unnecessary. 
Generally speaking, the zero-order data hold is found to be 
adequate for most systems found in practice. This being so, 
the data reconstruction process can easily be implemented 
by the digital computer since an output register of the com- 
puter can be utilized as a zero-order hold. 
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2. The effect of sampling. 

Sampling is the act of examining a continuous function 
at specified increments of an independent variable, usually 
time. Sampling processes can and do assume a wide variety 
of forms and can be separated roughly into three categories 

1. Linear, fixed pattern sampling 

2. Non-linear, signal dependent sampling 

3. Random sampling 

Included in the first category is the historically 
conventional concept of "instantaneous" sampling of fixed 
period T. By instantaneous it is meant that the duration, 
or pulse width, of the sample of the signal is very short 
compared to the shortest time constant of the system* Un- 
der this assumption (which will be employed in this thesis) 
the sampling process can be regarded as generating a se- 
quence of impulses. This viewpoint leads to the develop- 
ment of the Z-transform calculus and also to the use of 
difference equations as a mathematical tool. Extensive 
tables of Z-transforms may be found in the literature, and 
it is relatively easy to take the Z-transform of simple, 
low-order transfer functions. However, for a fifth order 
plant or higher, finding the Z-transform is quite laborious 
Appendix 2 describes a digital computer program written for 
transfer functions of Type 0 and Type 1 servo-mechanisms 
up to the tenth order. 

The Sampling Theorem by Shannon states that a band- 
limited signal of highest frequency F q can be completely 



5 



specified by sampling at a rate of 2F q . Recovery of the 
signal at this minimum theoretical sampling rate requires 
an ideal filter. In practice, signals are not ordinarily 
band-limited and filters do not have ideal attenuation 
characteristics. Even at high sampling rates, sampled- 
data systems usually exhibit distortion in the form of 
"ripple". One means of combating this ripple is to employ 
a multi-rate digital controller whose output occurs at a 
faster rate than the basic error-sampling rate of the over- 
all system. By applying control signals to the plant more 
frequently much of the ripple can be reduced and the plant 
can be made to respond more quickly to an input signal. 
Design judgment must be used, however, since the output 
signals from a multi-rate controller are frequently of 
greater magnitude than those from a single-rate device, 
and the possibility of plant saturation as well as the 
added complexity of the controller must be balanced against 
the advantages to be gained from a multi-rate system. 

Of crucial importance in a sampled-data system is the 

effect of sampling on the stability of the system. In an 

otherwise stable system, sampling reduces the degree of 

/ 2 / 

stability and can even lead to instability. ' There are 
various criteria for continuous systems that can be applied 
to sampled-data systems. Among these are the modified 
Routh-Hurwitz criterion, the Schur-Cohn criterion, and the 
Nyquist criterion. As an example of the determination of 
stability using the modified Routh-Hurwitz criterion we 

will use the third order Type 1 servo-mechanism shown in 

6 



Figure 2. The plant is closed-loop stable with the gain 
indicated. We desire to convert the continuous system to 
a sarnpled-data system / so we first insert sampling and data 
reconstruction as shown. The matter of digital compensation 
will be held in abeyance while we investigate the effect of 
sampling only. 



2-e£0-ORDE<? 



CONiTl NJ UOO S 



T 



— k> 
X' 



HOLD 


Plant 


-Ts 




go, OOO 


1 - <£ 




s 




S( 5 + SS-)(5 + 70) 



c(t) 



J 



Figure 2. A third order sampled-data system. 



Suppose that a sampling rate of 100 samples per sec- 
ond, or T .01 seconds, has been chosen a priori . Then 
the modified Routh-Hurwitz procedure gives a "yes or no" 
result in regard to stability. For this value of T the 
forward path pulse transfer function (or Z-transform) is: 

P(-B) .OlO&oo & Z +■ .033672. 2 + >0O6S-S'4 

Z 3 - 3. 31S38&Z Z + 7 2 -0. 3^67^/ 

For the unity feedback configuration shown the closed loop 

transfer function has a characteristic equation given by: 

P(z) + Q(z) = A(z) = A3 z^ + A2 z^ + Al z + AO 

where A3 = 1.000000 

A2 = -2.264786 
Al = 1.695799 

A0 = -0.380157 

The closed loop sampled-data system is stable if the roots 
of A(z) lie within the unit circle in the Z-plane. By 
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employing the Bilinear Transformation, Z = 1 + w , the in- 

1 - w 

terior of the Z-plane is mapped into the left half of the 

w-plane. Once the modified characteristics equation is 

3 2 

obtained, namely: A(w) = B3 w + B2 w + B1 w + BO then 

the Routh array may be formed and the criterion for a third 

order system may be applied. This criterion is: 

Bl*B2 - B0*B3 is a positive quantity if all the roots 
of A(w) are in the left half plane 

In this example, B3 = 5,340742 

B2 = 2.428516 
B1 = 0.179886 
BO = 0.050856 

so B1*B2 - 30*B3 = + 0.165247 and the system is stable for 
T equals .01 seconds. 

If it is desired to decrease the sampling rate as much 
as possible in order to time- share the computer or to allow 
for increased computation time between samples, then it be- 
comes necessary to know the limit of stability as the samp- 
ling period T is increased. Unfortunately, the character- 
istic equation is transcendental in T, and we cannot solve 
directly for T at the stability limit. Iterative procedures 
most readily delegated to a digital computer must be em- 
ployed. Appendix 3 shows the development of the character- 
istic equation and the Fortran program for the CDC 1604 
digital computer which calculates the Routh-Hurwitz criter- 
ion for this third order system for one hundred values of T 
ranging from ,001 seconds to 0.1 seconds. From the results 
shown the stability limit on T is evidently near .025 sec- 
onds, since this is the last value for which the R-H criter- 
ion is a positive quantity. 
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Figure 3 shows the response of this un-compen sated 
third order sampled data system to a unit step input for 
three values of T. For T of .001 seconds the response is 
indistinguishable to the eye from the results obtained in 
an analog computer simulation of the problem. This would 
be a good value to choose for the time increment in a 
solution of the plant differential equation by numerical 
methods (see Appendix 1). It can be seen that the sta- 
bility of the system is degraded somewhat when T is in- 
creased to .01 seconds, but the system is still stable. 
This is the sampling period to be employed later when we 
introduce digital compensation to this third order system. 
A sampling period of .05 seconds results in severe insta- 
bility, as we could expect from the Routh-Hurwitz deter- 
mination. Digital compensation is a must in this case, 
since stability is a prime requisite. Whether or not we 
compensate a stable system depends on how satisfactory the 
response is without compensation. 
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3. The digital process. 

As previously mentioned, a sampled-data system is 
characterized by the fact that in certain areas of the sys- 
tem signals appear as discrete quantities or fixed numbers 
in a sequence at regular intervals of time. The most ob- 
vious point at which this occurs is at the digital computer 
or controller. The input to the computer is a sequence of 
numbers and the computer output is also a sequence of num- 
bers. We can thus regard the computer as a "black box" as 
shown in Figure 4. 



ElOTt] 


Ducita l 


E2[pT 2 ] 






Computer 


P = l j • ' * j 



Figure 4. Generalized concept of a digital computer. 

The input sequence, El 0 T J • consists of N discrete 
values separated in time by T^ seconds. The output se- 
quence may contain a different number, M, of fixed quan- 
tities and the period, T 2 , of the output sequence need not 
be the same as the input period. 

For simplicity let us first take T, equal to T^. 

Then at a certain time, pT, the general form of linear 
computation that the digital computer can perform may be 
expressed ass 

e2 i>T] + bi-e2[(f>-0T] +•••+- b M _ ( • el [Cp-w+i)t] 

- clQ • el [pf] t ai‘ el [(P-OT] £ MH -ei[(p-N + «)T] 
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This can be written: 



A 1-1 



el OtJ = 2_ etf-K 

M- o 



M-i 



Z b J e2 p-j 

J=< 



or, by Z-transform calculus: 



^Z(s) 

£1(?) 



Z^kZ 
1 + Zbj z 1 



which can be considered the pulse transfer function of a 
digital computer. If all the b^ are zero then we have what 
is called a finite memory process since the output at pre- 
sent time, pT, is the weighted sum of the present input 
plus earlier inputs up to some finite limit, N. If not all 
of the b coefficients are zero then a weighted contribution 
of earlier outputs is also involved in forming the output 
at time pT. Since each earlier output was a combination of 
both input and output sequences, ad infinitum , we now have 
what is called an infinite memory process. Another way of 
showing this is to perform long division on the transfer 

function to obtain new weighting coefficients, c^., for the 



input sequence alone. Only in the case where a^Z con- 
tains 1 - bjZ”- 1 as a factor does long division fail to 
result in an infinite series, showing the need for an in- 
finite memory of previous inputs. When such a transfer 
function describes the over-all closed loop process, we 
have what is known as a "non-finite settling time" process 
and the systematic error (or difference between input and 
output) can never theoretically reach zero. 

In the system configuration employed in this thesis 
we will use the concept of the digital computer transfer 
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function in two ways. First is the over-all pulse transfer 
function of the closed-loop process,, defined ass 

/?(&) ~ //JPUT SSOUSTAJC-e 
Our attention will be devoted mostly to finite memory K(z) 

consisting of numerator polynomial only. Such a process 
has a finite settling time and can "forget" the effect of 
initial errors or acquisition transients in the input sig- 
nal. 



Second is the pulse transfer function of the digital 
controller in the forward path. Storage will be allocated 
for a finite number of values in both the input error se- 
quence and the output error sequence, and the controller 

-1 

transfer function will be a ratio of polynomials in Z 

defined as: _ , x _ Eout(z) _ processed error sequence 
' z ' Ein( z) observed error sequence 

In case the input sequence to the controller has a 
period, T^, different from the period, T 2 , of the output 
sequence we can still write a pulse transfer function that 
is amenable to interpretation using difference equations 
providing and T 2 are related by an integer. The pulse 
transfer function is then expressed in the "Z-domain" 
created by the faster sampler, which in this thesis will 
be at the output from the controller. In the next section 
is provided the mathematical basis for these assertions. 
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4. Response criteria. 

In view of the flexibility possible with active digi- 
tal compensation there is a very large number of possible 
over-all response functions which can be implemented. As 
a starting point, the simplest or "minimal " prototype re- 
sponse functions are convenient. Minimal prototype systems 
are approached with the point of view that they must be 
able to respond satisfactorily to some convenient test in- 
put such as a step, ramp, or constant acceleration, or all 
three. The requirements which are set for minimal prototype 
response functions are: 

a. The over-all response and the response of all ele- 
ments of the system must be physically realizable. 

b. The steady-state response to the test input must 
have zero systematic error. 

c. The transient response should be as fast as pos- 
sible and the settling time should be equal to a 

/4/ 

finite number of sampling intervals. 

The following table shows the test inputs and the required 
minimal prototype response functions, K(z): 



Input 




K(z) 


Step 


-1 

z 

^ -1 


-2 


Ramp 


2z 


- z 


Acceleration 


3z 


, -2 , -3 

- 3z + z 



These minimal responses are possible only for those plants 
whose pulse transfer functions, G(z), contain no zeros or 
poles outside the unit circle in the 2-plane. 
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The essence of digital compensation is that the con* 
stant coefficients of the pulse transfer function, D(z), 
for the controller are adjusted so that the desired over* 
all pulse transfer function of the process, K(z), is real- 
ized. Depending on the structure of the sampled-data feed- 
back system, the controller transfer function is derived 
from system functions in various ways. The configuration 
employed in this thesis is the typical one shown earlier in 
Figure 1, viz., in-line digital compensation of a unity 
feed-back, error-sampled control system. The ADC can be 
thought of as a single sampling switch operating with per- 
iod T. Data reconstruction will be accomplished by a zero- 
order hold, and the pulse transfer function of this DAC in 
combination with the continuous plant will be denoted by 
G(z). Since the actual input and output are continuous 
signals, fictitious samplers are shown producing R(z) and 
C(z) so that we may speak of the over-all pulse transfer 
function, K(z). For convenience. Figure 1 is re-drawn 
embodying these ideas as Figure 5 below. 
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Figure 5. Sampled-data system with in-line compensation. 



Taking the single-rate case first (where n-1, thus T/n-T), 
the over-all pulse transfer function is seen to be 
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K(z) 



C( z) D(z)«G(z) 

R(zT 1 + d(z)«g(z) 



which can be solved to give 



D(z) 



1 K(z) 

gTTT * i -, ”k(z7 



Depending on the poles and zeros of the plant which 

it is desired to compensate digitally* complete freedom of 

choice of K(z) may not exist. According to Franklin and 

Ragazzini (4) the restrictions on K(z) ares 

It is necessary that the specified over-all pulse 
transfer function, K(z), contain as its zeros all 
those zeros of the plant pulse transfer function, 
G(z) which lie on or outside the unit circle in 
the z-plane, and that l-K(z) contain as its zeros 
all those poles of the plant transfer function which 
lie on or outside the unit circle in the z=plane. 

Systems which are designed for minimum and finite 
settling time often do not give good performance when sub- 
jected to an input other than that for which they are de~ 
signed. In this sense, minimal systems may be regarded 
as highly "tuned". In addition, the severe shocks which 
result in the plant cause substantial ripple in the con- 
tinuous output even though the output is error-free at the 

sampling instants. These effects were noted in the lit- 
/7 / 

erature and led to the introduction of a term in the 
over-all response function known as the staleness factor. 
This factor leads to a softening of the response, with the 
result that the system can be expected to respond slightly 
slower but adequately to a number of test inputs. The 
choice of staleness factor can be arrived at by optimizing 

/g / 

procedures' ' or by observing the response to some most 
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likely form of test input. The staleness factor is intro- 
duced as the constant "c" in the following expression for 
the over-all prototype transfer function: 

tw_\ _ a»K( z)rainimal 
Z ' ~ —IN 

(l-cz-V 

where in practice, N = 1, c is a positive constant less 
than unity, and "a" is determined by applying the Final 
Value Theorem. 

In order to compensate the system so that the output 

is ripple-free after a finite settling time the following 

/ 4/ 

design rules must be applied;' 

a. All the rules for minimal prototype response 
systems apply. 

b. The open loop transfer function must be capable 
of generating an output that is the same as the 
input. 

c. The over-all pulse transfer function, K(z), must 
contain as its zeros all the zeros of the plant 
pulse transfer function, G(z), and not just the 
zeros of G(z) which lie outside the unit circle 
in the z-plane. 

Let us turn now to the case in which the output of 
the digital controller occurs more frequently than the in- 
put. The development that follows is due to Kranc (Refer- 
ence 9 ) . 

If we apply a discrete sequence, E(z), to a continuous 
plant, H(s), then it follows that the continuous output is; 

C(s) = H(s)-S(z) (1) 

For n equal unity we may take the pulse transform of both 
sides of equation (1) for a sample period of T. Since E(z) 
is already a function of Z it is invariant under the trans- 
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formation, so the result iss C(z) = H(z)»E(z) 

When the output is sampled "n" times more frequently 
than the input (n = 2,3,.. 0 ) let us adopt the notation that 
Z n denotes a pulse transform taken for a sampling period of 
T/n. Then we may defines 



°< Z n> 4 Z;{ c < s >} =£*{»<*> E < z >} 

H ‘ z n> S {*<=>} 

Once again E(z) is invariant (also see Reference 4, page 88) 
so C(Z n ) ~ H(Z n )*E(z) but we have a mixture of z and Z^ 
transforms, and the meaning is unclear,, If we now define 



then it can be seen that Z = zJJ and E(Z) - E(zjj) so in the 
Z n "domain" we can write C(Z n ) = H(Z n ) E(Z^) C (2) 

Returning now to the closed loop system of Figure 5 
we can see the following relationships? 

*<*>’ %8r, = ^ 5 ?) 

E1(h) = R(z) - CCz) (3) 



C (z) - El (?) (Zn)‘(? ^2-n) j- (4) 

C (2*) = E1(2)-P(2v.)* G-CZy,) (5) 

D/g ) 

Using (3) and (4): El(2) = — r (6) 



Substituting 
(6) into (5): 



r(z) • 2LEz±1 : (?) 

1 + &(*■»)} 
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Since Kdy,)- there results: K(2-») = — — : 

xt**) i + Zfpt&ycG*)} 






which can be re-written as: 

KteJ/C {&(*■”)• £<■*”)} -Ote»). £•(*■») 4 K(2»)=0 

Equation (3) is a functional equation of the form 

Z{F(a„)} - F(2„) = - K (2-k ) 



(8) 



It can be shown 



/9/ 



that the solution is 



F(*0 = ' - D(2v,)' CrC. Z ^) 

so finally, the transfer function of the multi-rate con- 
troller is: 



D / 2 n - — ! _J <C g A _ 



(9) 



The Z n ~ transform of the plant is simply obtained by 
substituting Z n for Z and T/n for T in G( z). The Z - trans- 
form of the overall process is not so simply obtained. We 
must start by imposing steady-state criteria and the re- 
straint of physical realizability. The following develop- 
ment is due to Franklin and Ragazzini (Reference 4). 

For inputs whose Laplace transform is R(s) - — we 

S K 

can write P(2-«) 



RC^WttV k(^) = -E^ = 

( \ - Zy\ ) (\ ^ 



C(2v) 

(l-2^)K K(2„) 

Thus the error sequence, R(Z n ) - C(Z n ), can be expressed as 

uo) 

For a stable process the Final Value Theorem states that 

El (<*>) = , (i-C)El(^) UD 

Z-W + l 

In order for equation (10) to meet the requirements of 
equation (11) we must have 

K(Z*) = U2> 



(j-z; 1 )* 
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where f(Z ) is an arbitrary function of Z whose coeffic- 
n n 

ients are to be determined. 

Substituting (12) into (10) results ins 



s. PO*) - 

E1(2 " ) - (1-2,-r 

and for zero steady state error the numerator of 2l(z) 
must contain as a factor the denominator of the inputs 



(13) 



R(z; n ), i.e., 



0-^') K -T(2^ = Q(2„) (14) 



or equivalently: 



=■ O -for WisO,i ...K 
2 -m= 1 U5) 



The restraints imposed by equation (15) are applied in order 
to determine the coefficients of the arbitrary function,, 
f(Z n ). Once f( s n ) is known, equation (12) gives K(Z i ). 

The only unknown now remaining in the fundamental 
design equation for the digital controller , equation (9), 
is K(zj^). This is obtained from K(Z n ) by merely taking the 
j powered terms of K(Z n ). . ..j = n, 2n, 3n, ... 

In very special circumstances, K(Z^) may be found in 
a simpler way by using so-called "single-rate 1 ” techniques 
described in Reference 4, These special circumstances ares 

a. The plant must be open-loop stable. 

b. K(Z n ) I s limited to a finite settling time type 



c. 



(i.e., numerator polynomial only). 

The speed-up factor, "n‘ ! , in the controller out- 
put must be selected so that all transients at 
output sampling instants are over in one period 

of the input sampler, T. This is equivalent to 
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requiring that "n" equal the number of terms in 
the numerator of the plant # s pulse transfer func- 
tion, G(z). 



If these circumstances obtain, and if the design goal is 
to achieve ripple-free response to a step input, then for 
a plant transfer function of 



The main limitation of this design procedure using 
single-rate techniques is that the speed-up factor, "n", 
must be the same as the number of numerator terms in G(z), 
which may be out of the question for a high order plant* 

For the fifth order plant to be discussed later a double- 
rate controller gave improved performance, but this con- 
troller could not be designed by the above procedure* It 
is also worth noting that the number of terms in D(Z n ) can- 
not be less than "n". 

Up to this point the discussion of response criteria 
has been predicated upon the input being a deterministic, 
noise-free signal, but it also may be desirable to judge 





/ 4 / 

it can be shown' that the digital controller transfer 



function is given by: 
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a process on the basis of its noise rejection capabilities. 
The following remarks concerning noise rejection are due to 
Monroe (Reference 1). 

In order to consider an input signal which is contam- 
inated by noise, certain simplifying assumptions are neces- 
sary. First, assume that the system in the absence of noise 
would perform the ideal operation on the input, ie. , the 
systematic error is zero. Second, assume the signal and 
noise are uncorrelated. Last, assume the noise is "white" 
noise over at least one sampling period. Then it can be 
shown that the mean squared error of the process is: 



The "variance reduction factor" is defined as the 
ratio of the mean squared error of the output function to 
the mean squared value of the input noise and is seen to be 
the sum of squares of the coefficients of the digital pro- 
cess. Clearly, the smaller this factor, the better the 
noise rejection capabilities of the digital program will be, 
and in general it is desired that this factor shall be less 
than unity. In passing, note that the minimal prototype 
responses to ramp and acceleration inputs have variance 
reduction factors of five and nineteen, respectively. 

The digital process can be constrained to have no 




*■/ 



% 

where is the mean squared value of the input noise 
CL ^ is a coefficient of the digital process, K(z) 

Kte) = <2k 
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systematic error when the input is taken as a polynomial 
of a given order or lower. For example, if linear con- 
straints are imposed, then no systematic error exists for 
ramp or step inputs. If, subject to linear constraints, 
the variance reduction factor is minimized, the following 
expression is obtained for the coefficients of a digital 
process employing N samples of past history and predicting 
forward Od sampling intervals. 



The requirements that K(z) be able to predict for- 
ward byodT arises due to the delay inherent in the computing 
process, D(z). The discrete error value at sampling time nT 
equals the input signal at time nT minus the output result- 
ing from a control pulse applied to the plant one sample 
period earlier, ie., at time (n-l)T e In order to have zero 
error, the output must equal the input, so K(z) must predict 
forward one sample period, thus ©< = 1. 

The second formula above can be solved for N, given 
any desired value of variance reduction factor. In the case 
of unity variance reduction the result shows that six samples 
are required, and the process pulse transfer function is: 



; -h & << __ (a/ - / -a a od j 

A/ ( A/ + / ) A/ ( A/ + > )( aJ - / ) 



and the variance reduction factor is: 




#= c 




— 
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It was noted earlier that the presence of a plant zero 



on or outside the unit circle in the z-plane prevents the 
realization of minimal prototype response. Similarly, such 
a zero will prevent realization of the unity variance re- 
sponse. Indeed, when a zero must be accounted for, it be- 
comes quite laborious to determine any set of coefficients, 
a, , that will achieve noise reduction. According to Monroe 



(Reference 1) it is possible to design K(z) with noise re- 
ducing properties for a plant with a single zero outside 
the unit circle by the following procedures 

let b = the value of the zero to be accounted for. 



The pulse transfer function of the process using N data 
points is y_ t 

k(*)= O' 

For linear constraints on the process (ie., zero systematic 
error for ramp or step inputs) then the matrix equation be- 



then employed in the recursion formula shown in order to 



= unity for prediction forward on one sample 



period. 




These values are 
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It was noted earlier that the presence of a plant zero 
on or outside the unit circle in the z-plane prevents the 
realization of minimal prototype response. Similarly,, such 
a zero will prevent realization of the unity variance re- 
sponse. Indeed, when a zero must be accounted for, it be- 
comes quite laborious to determine any set of coefficients, 
a^, that will achieve noise reduction. According to Monroe 
(Reference 1) it is possible to design K(z) with noise re- 
ducing properties for a plant with a single zero outside 
the unit circle by the following procedures 

let b = the value of the zero to be accounted for. 



The pulse transfer function of the process using N data 
points i « 



For linear constraints on the process (ie., zero systematic 
error for ramp or step inputs) then the matrix equation be- 



then employed in the recursion formula shown in order to 



oC = unity for prediction forward on one sample 



period. 





These values are 
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The recursion formula is: 



0 
f M 



H +4 




(/-b+b 1 )* 

i rxi-T* 




Sven with a noisy input signal it is still desirable 
to exercise some control over the transient behaviour of 
the system. It is possible, again according to Monroe, to 
design the over-all process according to a composite cri- 
terion wherein preference can be expressed for step tran- 
sient response (in the mean squared error sense) over noise 
rejection by means of an arbitrary parameter , ft , ranging 
between zero and one. Subject to linear constraints and 
for no plant zero outside the unit circle, the process co- 
efficients, a^, are calculated from the following formulae: 
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(^recursion fe/cL*k/ ‘oyiship ) 



#i+X = 



tJ(i-fl) 

tf . 0 — ^dl. ~ ^ - Ao *4 — A / ^ (^ £~£a.r't!-y\cj 



A », -A,l>* 



A. 4 - A ^ 



*Fc rmu lae ) 






4 = Ahi 

X,-' 



M, =. (»-‘)A?'- NX? + X 



A/ 






_ (N-/) X^- A/ A 7 ■< L X/ 



a; 



J. 

2 



^x 0 -/j x f a, -/; x 

A ° = x f 2 + f % 7 ^} ] 

X[= z( z } ~ i ^Hth ) ] i 



Figure 6 shows the Fortran coding of these equations 
and the results for a ten point process with beta of 0.75. 
The program may be used for another design merely by chang- 
ing the two statements giving the values of the parameters 
beta and AN. 

For the design of a process to a composite criterion 
when the plant has a zero outside the unit circle see Ref- 
erence 1. 
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JOB STAPP3 

PROGRAM STAPP3 

COMPUTES COEFFICIENTS OF THE PROCESS KIZI TO COMPOSITE CRITERIA 

B IS A PARAMETER SHOWING THE PROPORTION OF EFFORT OlRECTEO AT i. 
INTEGRAL SQUARED ERROR TRANSIENT RESPONSE BEHAVIOUR 

!n°S 7 ?H6 NUMBER OF DATA POINTS TO BE USED IN THe PROCESS K U) 
AN* 10.0 

C*B/ ( AN» ( 1 .0-0) ) ; 
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Figure 6. Progrcun for computing composite coefficients 
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5. Program requirements. 

Before discussing the realization of the digital com- 
pensator, D(z), by programming of a digital computer, let 
it be noted that the pulse transfer function of the compen- 
sator can be realized in at least two other ways if a digital 

/ 2 / 

computer is not available or economically feasible. First 

is the realization by a delay-line network, sometimes refer- 
red to as a pulsed-data processing unit, through the use of 
delay elements, potentiometers, and summing amplifiers. The 
pulse transfer function can also be realized by a pulsed- 
data network consisting of resistors and capacitors. In- 
ductors are avoided, because in the usual frequency range of 
digital and sampled-data control systems the sizes and weights 
of the required inductors would be prohibitive and the losses 
would be excessive. Pulsed data RC networks have the advan- 
tages of simplicity in structure and economy in implementation. 

When a digital computer is available, realization is 
relatively simple by programming. The program can generally 
be carried out by means of one of three different methods: 

a. Direct programming 

b. Iterative, or cascade, programming 

c. Parallel programming 

Let the pulse transfer function (or Z-transform) of 
the digital controller be 
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The first terra in the denominator of D(z) must be a non- 
zero constant for the sake of physical realizability. In 
the direct method of programming it is necessary that this 
term be unity, since the direct method employs the differ- 
ence equation equivalent of the above pulse transfer func- 
tion. Writing the difference equation we get 
N-l M-l 

EoGT [’ nT l - Em ~ 2i Ecur 

In the digital computer this method requires 2M+2N-1 
cells of storage for data and constants, M+N-l multiplica- 
tions, M+N-2 additions and/or subtractions, and M+N-2 data 
transfers. 

The cascade method is based on the fact that D(z) is 
factorable into its poles and zeros such that a cascade of 
first order transfer function can be formed. In the computer 
this method requires 2M+N+3 cells of total storage, M+N+l 
multiplications, M+N additions and/or subtractions, and 
M+2 data transfers. 

The parallel method of programming is based on the 
partial fraction expansion of D(z) and requires 3M+2 cells 
of storage, 2 M multiplications, 2M-1 additions and/or sub- 
tractions, and M+2 data transf ers.^^ 

Before a choice of methods can be made on the basis 
of storage requirements or computation time, D(z) must be 
known, since any advantage of one method over the others 
depends primarily on the number of terms in the numerator 
of D(z). Generally, other considerations prevail. For 
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example, in the cascade method each stored constant is either 
a pole or a zero of D(z), and thus this method lends itself 
readily to experimental manipulation of the pulse transfer 
function in a "cut and try" sense. If storage space is at 
a premium, as in a special purpose digital computer, the 
direct method enjoys the advantage that data transfers can 
be effected by a single command that processes all the data 
in a shift register . In the simulation program described 
in Appendix 1 the direct method of programming the digital 
controller is used. 
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6. Implementation of a double-rate digital controller by 
direct programming methods. 

Direct programming takes the difference equation 
equivalent of the pulse transfer function and interprets 
z" 1 to be a delay of one computation period. Thus, in the 
single-rate controller both input and output discrete se- 
quences are equally spaced in time, and the programming is 
straight-forward. The computer program for a double-rate 
controller will be described in terms of the various arith- 
metical and shifting operations that must take place and 
not in terms of any particular programming language. The 
method will be shown using an example, and can be easily 
extended to the case of an n-rate controller. 

Suppose that the pulse transfer function of the double- 
rate controller is: 

0(2 Z ) - 4 AT*'?* + ^3 +• + A S ' • + AA 1 

Z ~ 1 + +- + 6 + Bte'Zl'*' 

The output from the controller will be a sequence, SoutCZ^.), 
of discrete values separated in time by T/2 seconds. The 
input sequence is spaced at intervals of T as determined by 
the error sampling rate. When we write the difference equa- 
tion to find a new value for Eout by cross-multiplying the 
transfer function and treating as a ^ e l a y °f T/2 seconds 

we get: 

eout = Al- + AZ • Ein 

+ A 4 1 £ 'I s1 + &£• Eim + Ab' Eito 
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i 



BZ*Eout 

+ BS“‘ Eout[(n- 4) t / ji ’] 



-h 83 • Eoot [_ (H-z) r /z~j +■ BA Eout [^M-3)1? a ”] 
4- 06 i Eoot jY N j 



Since the discrete sequence, Ein(z), is at intervals of T 
seconds only half of the terms involving Ein above will be 
non-zero. Which half depends on whether N is even or odd. 
When N is even (ie., an output value is generated at the 
"same" time the error is sampled) then the terms involving 
Al, A3, and A5 are non-zero. When N is odd (ie. , an output 
value is produced by the controller between input samples) 
then the only terms with non-zero values of Ein are those 
involving A2, A4 , and A6. This amounts to alternately 
multiplying the input sequence by the even numbered coef- 
ficients of the numerator of D^i^), then the odd numbered 
ones, ...etc. 

In a given period of time the output sequence from a 
double-rate digital controller will consist of twice as 
many discrete values as the input sequence. For the sake 
of illustration let us number both sequences backward in 
time in the same way, thus: The most recent value is num- 

ber one, the next most recent value is number two, etc. 

The total number of values in each sequence that must be 
stored in the computer memory is determined by the number 
of terms in the controller transfer function. In the ex- 
ample of this section we require output values numbers one 
through six and input values one, two and three. Each out- 
put value is employed only once in each numbered position 
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and passes from storage after six intervals of T/2. Each 
input value is employed twice in each numbered position ^ 
once as a multiplier of Aj and T/2 seconds later as a mul- 
tiplier of A.,,. Since the numerator of D(Z_) contains 

J TJL £ 

six terms, each input value also passes from storage after 
six intervals of T/2. A simple-minded pictorial represen- 
tation may serve to clarify these ideas. Let the present 
time be the 16th instant of the closing of the output samp- 
ling switch of the double-rate controller,, Neglecting 
computation delay (which is usually at least an order of 
magnitude less than the sampling period) the present time 
is also the eighth closing of the input sampling switch. 

The input and output sequences may be thought of as exist- 
ing in storage as shown in Figure 7. 

Input sequence , Ein 



observation 

©^ 



© 



© 



ST 



present 

time 



7T 



GT 



time 



Output sequence , Sout 

® © 



© 



to be computed 

/ t t 

16T/2 15T/2 14 T/2 13T/2 12T/2 11T/2 



time 



(older values) 

Figure 7, Input and output sequences for EK^) 



In the difference equation the terms involving the input 
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sequence would be (since in this example, N = 16): 
Al*Ein(8T) or Al*Ein(l) 

A3*Ein(7T) or A3*Ein(2) 

A5‘Ein(6T) or A5*Ein(3) 

The value of Ein(l) has just been measured by the input 
sampler and the controller produces Eout(l). 

Now let us move on in time by T/2 seconds. It is 
necessary to re-number the output sequence since a new 
Eout(l) will be produced. The previous Eout(l) becomes the 
current Eout(2), etc., with the previous Eout(6) being dis- 
carded. We are now between input samples and no new input 
information is to be received. The previous input values 
have become T/2 seconds older, however, and the sequences 
in storage might be pictured as shown in Figure 8 e 

Input sequence , Ein 






(D 



8T 

to be computed 

d/ ® 



17 T/2 IS T/2 

present 

time 






7T 

Output sequence , Eout 

© 0 



6T time 






15T/2 14T/2 ~3T/2 “ T2T/2 

(older values) 



time 



Figure 8. Sequences between input samples 



Since the difference equation is in terms of T/2 the non-zero 
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terms involving the input sequence have changed and are 
now (for N = 17): 

A2»Ein(8T) or A2»3in(l) 

A4* Ein ( 7T) or A4*Sin(2) 

A6*Ein(6T) or A6*Ein(3) 

The controller computes a new output value using all 
of the denominator coefficients of its transfer function 
but only the even-numbered numerator coefficients. As 
before, the output sequence must be shifted since we are 
only half a period away from receiving a new observed value 
for Ein(l). 

For an example of the implementation of these ideas 
using the CDC Fortran programming language see Appendix 1. 
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7. Applications. 

The results to be described in this section were 
obtained using the digital computer program simulation of 
a control system described in Appendix 1. The response 
curves were obtained through use of a Graphplot computer 
routine for the CDC 160 digital computer devised by Lt. 
Robert L. Hogg of the U. S, Naval Postgraduate School , 
Monterey, California. The program in Appendix 1 also pro- 
duces numerical results on a printer, but these are not in- 
cluded in this work. 

Initial investigation was made of a simple plant 
described by a second order differential equation whose 
pulse transfer function had no zero outside the unit circle 
in the z-plane. Next, a third order plant with a pulse 
transfer function containing a zero outside the unit cir- 
cle was investigated. Both of these plants were theoretical 
and in neither case was an actual physical counter-part at 
hand. 

At the U. S. Navy Electronics Laboratory, San Diego, 
California, a group headed by John B. Slaughter is engaged 
in studies of Digital Sampled-Data Feedback Systems for 
Shipboard Control Applications (NEL problem N4-17). The 
author was privileged to spend ten weeks with that group 
early in 1963, during which time a twin 3"/50 gun mount 
was in the process of being installed to serve as the plant 
in an experimental control system. A USQ-20 general purpose 
digital computer is to be employed as the controller. The 
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manufacturers of the gun mount supplied a fifth order 
transfer function for the gun in train. The majority of 
the simulation studies to be described in this thesis were 
conducted using this transfer function, and it is felt that 
the results obtained are applicable not only to this par- 
ticular obsolescent plant but also to present day ship- 
board weapons, especially missile launchers. 

The simulation program has the capability of testing 
a system for input signals of either unit step, ramp, or 
acceleration. In addition, additive Gaussian noise of 
selected mean and variance may be included m the input. 

The first action of the program is to produce the response 
obtained in the case of the continuous, un-compensated sys- 
tem for purposes of comparison. In the event this response 
is unstable no graph is made of it. Instead, the test input 
is reproduced. The next action of the program is to simulate 
the same plant incorporated in a sampled-data system with a 
single-rate digital controller and produce the response for 
this case. The program is further capable of simulating a 
system using a double-rate digital controller to produce a 
third response curve, although this is not done in every ex- 
periment. For convenience, the program will also compute 
the pulse transfer function of the controller when requested 
to do so, or the controller coefficients may be provided by 
the user. The Z-transforms employed in this work were com- 
puted by the program described in Appendix 2 and were veri- 
fied by hand for the second and third order systems. 
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The second order plant has a Laplace transfer func- 
tion given by i / /^\ _ // 41 

HCS ~ S(s-hZS) 

The pulse-transfer function of the plant and zero-order hold 
combination for a sampling period of T = .01 seconds is 

(l--2~')(/-6-~71 2*0! Z") ~ /- /. 77220/ H'V 0.77$qoi 2~‘ 

There are no zeros of G(z) on or outside the unit circle in 

the z-plane, so minimal prototype compensation can be a~ 

chieved. As noted earlier , the minimal prototype over-all 

pulse transfer function for a step input is K(z) = z“^ , so 

I ^( 2 ) __ /- /,-77 ggo/e‘V .JJJJCjl g~' 

' <?(2) I - « (2) ~ . as’222~‘ +- *64 86 2 2~‘ 

2 _ /£73f3fry-/# 7^40/9 2"' 
• Og'L%^ r (l + . c il.G>4SSZ~') /.O -h 0.9204 S’-TZ"' 

(note that l-K(z) contains the pole of G(z) lying on the 
unit circle. ) 

Figure 9 shows the plant behaviour in a continuous, un- 
compensated system and the behaviour in a sampled-data 
system using the digital controller above. Note that there 
is zero systematic error at the first sampling instant after 
the input was applied, and also at all succeeding sample 
instants. For a K(z) containing "n" terms there will be 
zero error for a design input after "n" sampling instants. 

The plant behaviour between samples is another matter, and 
in this example probably would not be considered satisfac- 
tory, unless the only desire was to increase the speed of 
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response with no regard for over-shoot and inter-sample 
ripple. 



Simulated Gaussian noise with a standard deviation 
of 0.204 was added to the input and the experiment was re- 
peated. Figure 10(a) shows the results. According to Mon- 
roe the process should have unity variance reduction factor 
since the sum of the squares of the process coefficients is 
one. Judging from Figure 10(b), which shows the noisy test 
input, this is approximately correct, since no noise reduc- 
tion is apparent. 

In an effort to improve the transient behaviour for a 
step input a staleness factor of 0.5 was selected. The 



over-all process was designed from 



K(*) - 



a- YU rrtt 'yntL / 



/- O.S 2 "' 

and the Final Value Theorem showed that for zero steady- 



state error, a = l- c = C.5 so that K(s:) = 



O. 



-/ 



/ - O. S’. 



D(2) = 



, OS-Z-' 

(/-■Z 'J(/-0.77S*0/S J — 



. OS2X E~'+ . 04 2 1 S' 



/ 



/- O.^TZ 



_ O-S'- 0. 3 29400^2 
S- . 04202.-' 

The step response using this digital controller is shown 
in Figure 11. This is an example of a non-finite settling 
time process, but for practical purposes zero error is 
achieved at about the seventh sampling instant and there- 
after, and the inter-sample ripple and over-shoot is much 



<2 40W7 — 7 374‘aa 
/*Q a O. 930 4 S' s' Z 
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reduced compared to the minimal design. 

This same staleness factor design was subjected to 
a noisy step input with the results shown in Figure 12. 

Since K(z) can be expressed as an infinite series by long 
division, thus: 0.5z -1 + . 25z~^ + .125z~ 3 + . 0625z~ 4 . . . . 

the variance reduction factor is approximately one-third. 

It is not apparent from the simulation results that the 
process indeed had this noise reducing property. 

The six point finite memory process for unity variance 
reduction (see page 23 above) was then subjected to a step 
input with and without noise. Figures 13 and 14 show the 
results. This process and the following one depend on the 
pulse transfer function of the plant not having a zero out- 
side the unit circle. It can be seen in Figure 13 that the 
transient behaviour of a design whose sole criterion is 
variance reduction is rather poor. For this reason the next 
process considered v/as a composite one in which three-quarters 
of the design effort was towards transient control and only 
one-quarter towards noise reduction. The program on page 
27 was used to compute the process coefficients, as shown 
on the same page. This gave the most complex controller so 
far with 11 numerator terms and 11 denominator terms. The 
results of the simulation experiment with this controller 
are shown in Figures 15 and 16. The transient behaviour 
displayed in Figure 15 is a noticeable improvement over that 
of Figure 13 where the only aim was noise reduction. 

At this point a certain disenchantment with the notion 
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of noise reduction begins to set in. The plant is a simple,, 
well-behaved one, and in the four experiments so far using 
"clean’ 1 input signals both the continuous system responses 
and the sampled-data system responses can be verified ana- 
lytically with some effort. When the unit step input has 
"white" noise added to it as shown in Figure 10(b) two ten- 
tative conclusions suggest themselves. First, the contin- 
uous system with no compensation is relatively unaffected 
by the presence of the noise. A close inspection of Figure 
10(a) reveals that the peak overshoot still occurs at about 
one-tenth of a second and its value is increased only slight- 
ly (from 29% to 37%). A steady state error equal to the mean 
squared value of the noise can be inferred from Figures 12 
and 16 in which the time scale has been doubled for a longer 
look at the plant 0 s response. Over-all , however, the low- 
pass character of the plant effectively prevents any drastic 
reaction to the presence of the noise. Not so in the case of 
the sampled-data system, so our second thought at this point 
is that the insertion of the digital controller in the for- 
ward path disrupts the low-pass filtering action of the plant. 
The noise in the error signal is the same as in the input, and 
the noise power in each discrete error sample is integrated 
over the next sampling period and applied to the plant. This 
integrated noise power causes some fairly drastic gyrations 
in the plant regardless of the design criteria. A brief com- 
parison of the four noisy experiments up to this point provides 
small reason for employing a variance reduction design in 
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preference to a noise-free criterion. It appears that the 
effort could better be expended in pre-filtering or other 
measures to eliminate noise prior to the input to a sampled- 
data system. 

The next four simulation studies were conducted using 
a third order plant whose transfer function is given by 





s( s s- as-J(s+ 70 ) 



For T = .01 seconds the z-transform of the plant-and-hold is 






. O/O &O0-2 . 0336 72. -t- - O O <Z> g Z 3 



/.O - 3. 37S32C>-2~‘+ / 6i^72" X -^3?67^/2 3 
or in factored form 

a. 967 37S~ 2~ ' ) 

(/ - 0.7 7226/Z~'X/- B~') 

The presence of a zero of G(z) outside the unit circle is 
noted. This zero will prevent us from realizing a "minimal" 
prototype response. In accordance with the rules given on 
page 16 the design of a process giving prototype response 
to a step input proceeds as follows: 

- (/ + 3'90>727<s- 2~') 

!-/<(*) — (/- z-)(/+ b x z~") 

Solving for the two unknowns, and b^ # results in 

.3&06Z ^*.147 93* and thus 

K( 2) = O. 2S20G, HZ."' ~h O. 747 93 2 Z~^ 

= ft + 3. 70 7 3 7S“Z ~')0‘ 33-30622-' 

Solving for the digital filter coefficients by hand, we write 
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tf>Z^(/-,778$0/B')(/-.496S?S'2') JS20&2^ ' 

^ = ('.0/at)(/+. 309327j/Xte^ J¥ * : Kr7Z^C^ 

and after cancellations of the factors shown, we get 

9(2) = 23.719934 - 3Q.3^7f?3£~y 

/.O i- 0, 9S73&-7 E~' 9 0. /sT^ 2 2“ a 

The un-compensated response and the response of the sampled- 
data system using the D(z) above are shown in Figure 17. 

The non-minimal nature of the digitally-controlled response 
is characterized by the fact that it took two sampling per- 
iods (.02 seconds) for the output to reach zero error at 
sampling instants. There is a slight amount of inter- 
sample ripple, but the peak overshoot is greatly reduced 
from what it was for the continuous system. The improved 
damping of this third order sampled-data system compared to 
the second order system considered previously is attributed 
to additional filtering action of the pole at -70. Since 
this design was for a step input there will not be zero sys- 
tematic error for higher order inputs, such as a ramp. This 
can be observed in Figure 18 where the curve labelled D1 is 
the response of the digitally controlled system to a ramp 
input. (Note: in this figure, and all others for which the 

test input is a ramp, a 45 degree sLope from the origin is 
the input signal). The controlled response settles down to 
a steady state error of about .017 units after two sample 
periods. Whether or not this is an intolerable error in re- 
sponse to a velocity input depends on the system. For ship- 
board fire-control application such as tracking a target, a 
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design for zero error to a ramp is usually preferable. 

More about this later. 

An additional refinement of the step response of this 
third order system can be achieved at the cost of one extra 
period, T, in reaching steady state by designing for ripple- 
free response as discussed on page 17. For this case we get 

tffr) = (/ r O- 3093 27 Z~’)(/ + ^ 

- (/- *~'X/ y- + b& 2~ Z ) 

and by equating coefficients of like powers of z' k (k=l,2,3) 
we obtain the below equations to solve for al, bl, and b2. 

- / - b ± 

^ • / ~7 Cs Q 2 . ^ 2 ^ ~ ^ 

0. 62//30 2 O.JL = 

Solving these, we get al = .2084318 so that the process is 
K(z) = . 2084318z -1 + ,662105z“ 2 + . 12946 34z~ 3 

These process coefficients as well as the pulse trans- 
fer function of the plant were provided to the simulator 
program which then calculated the digital controller and 
produced the ripple-free response shown in Figure 19. 

Figure 20 shows the same system with a velocity input. The 
steady state ramp error has increased slightly, from .017 to 
.020 due to the additional term in the process z-transform, 
and the controller had a total of six terms. In the absence 
of input noise we have demonstrated satisfactory digital con- 
trol of this lightly damped plant when the input is a unit 
step. 
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Turning now to the twin 3"/50 gun mount*, the transfer 
function of the forward path train drive is given by 

/ / / g \ _ /• ^ O x 

5“ y- // Jfs* v- y- + * 3 ) 

Most of this very large forward gam is supplied by the 

amplidyne drive units and results in a quickly responsive 

system when suitable compensation is added. In actual use 

the train drive function was elaborately compensated by 

various continuous networks in a multitude of feedback 

paths in order to obtain stable operation. All of this 

continuous compensation has been removed from the gun by 

the group at NEL* and the necessary control will be achieved 

by digital methods. Figures 21 and 22 were provided by 

Donald Lackowski* an engineer in that NEL group* and show 

the root locus and Routh criterion analysis for the forward 

transfer function with unity feedback. It cart be seen that 

the stability limit on the gain* K* is roughly three orders 

of magnitude lower than the actual value present* and thus 

the plant is highly unstable without compensation. This is 

borne out by the simulation studies to be described next. 

In every case the continuous response went rapidly off scale 

and could not be shown in the figures to follow. As previ= 

ously mentioned* in this event the simulator program re 

produced the test input in lieu of the unstable response. 

The computer program described in Appendix 2 was ex- 
ceedingly helpful at this point since finding the z-trans- 
form of this fifth order equation for just a single value 
of T is quite laborious by hand. With the computer a 
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number of transforms was quickly obtained for a spectrum 
of sampling rates ranging from one hundred samples per 
second to as slow as one sample every four seconds* For 
the faster rates the pulse transfer functions themselves 
were of the fifth order (ie.„ five numerator terms and six 
denominator terms). As the rate was decreased it was ob- 
served that the higher ordered terms in both numerator and 
denominator of G(z) became progressively smaller* For T 
of one-half second or larger a good approximation of G(z) 
can be made using only a second order function* In addi- 
tion, the zero of G(z) which lay outside the unit circle 
for smaller values of T had moved inside the circle for 
T = 0.5 seconds. The computer calculated pulse transfer 
function of the plant and hold combination for T of one- 
half second is given below: 



Numerator of G(z) 


Denominator of 


1336.33153319Z" 1 


1.00000000 


897. 96714053z“ 2 


-1. 004300 37 z" 1 


6.02702748z~ 3 


*00430125z" 2 


. 00049 220z" 4 


-.0000Q087z~ 3 


.00000000z“ 5 


. 00000000 z** 4 
.00G000Q0z“ 5 



The computer program also factored the numerator and gave 
the following real values for the zeros of G(z)? 
-.00000001, -.00008267, -.00669675, -.66518491, 0*0 
It was decided to investigate the advantages of a multi- 
rate sampled data system, the analysis and design of which 
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was discussed beginning on page 17. The simulator pro- 
gram has the ability not only to produce results for the 
continuous system and a single-rate digitally controlled 
system, but also, upon request, to simulate a system using 
a digital controller whose output to the plant occurs 
twice as frequently as the input error samples are sup- 
plied to the controller. This is called a doub le- rate 
controller. The output rate was chosen to be two samples 
per second in order to take advantage of the second order 
approximation to G(z) mentioned above. It was desired 
that the over-all process have zero steady state error to 
a velocity (ramp) input. Following the method of page 19 
the design for prototype ramp response proceeds as follows s 
The fundamental design equation for the double-rate con- 



Applying steady state and realizability constraints, we get 



troller is 




For T/2 = 0.5 seconds, G(Z y ) ~ g — 

^ / —• if ^2 





and a look ahead at the constraint equation, 








Hj'l 
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shows the need for only two free constants, a, and a-* so 

0(2*) = " Pfei. )* F(^x) -T*£ a (Q±2z' + 

= o.rzj'-z^&zr+^z;*; = o.&;‘- a^-*. a 3 z\ 



■for m = 0 : QCH^) 



2, = 1 



- 0,S 



C !j_ — <0^ — o 



■for ~m = 1 t dQ 

a?;' 



= o.s- - ?a±- 4 a*. = o 



Za*l 



solving these equations we gets 



a — 1 . 5 , a 



- 1.0 



so that K(z a )= (i^J ) (l.rz-'-2r)= (i + z^i-rej'-z; 1 

= /. re/V j.az' 3 - o.sz' 3 ~ /.a s ~ 4 

To find Ktei) we take the even -powered terms with the re- 

suit K(z*)= 3. 0 - /. o z~ 4 = Kfei) 

which can be recognized as the minimal prototype ramp 
response function shown on page 14 expressed in the 
domain. 

Substituting now in the fundamental design equation, 
the pulse transfer function of the double-rate controller 

n , x (/' 2*) (/• * 2-° ** a -a^ga 3 - 

^ 2 / (/534>b;'+ mzl*j ( / - 3 . O 2 " 3 -h /.0 4) 

= / - gj 7 

/336 - 43 2 Z^ - 2 “* 

We must normalize this result with respect to the leading 
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denominator coefficient, since direct programming requires 
a value of unity in this position. The final result is 



DteJ = 



.60//a27S-4 - . 0007 4 $ S O 
/.0 - 0 . 327 2443 ! — 0 . 672/^69 Z " x 



Before proceeding with the simulation run we must 
provide for the single-rate action with T = 1.0 seconds. 
This is simply done by providing the simulator program with 
the ramp minimal prototype function and G(z) coefficients. 
For this sampling rate an even better second order approx- 
imation is available for G(z), namely? 

, , , _ 3S-7& 3~' + 924 g~ a 

6rt * J_ , _ _-/ 



The responses of the two systems for a noise-free 
velocity input can be seen in Figure 23. The single-rare 
controller gave zero error at the second sampling instant 
and only moderate overshoot and inter-sample ripple. The 
double-rate system was 0. 5 seconds faster but over-shoot 
was worse due to the larger control forces involved. The 
ripple in the double-rate case was more persistent, also. 

These same systems were then subjected to step inputs 
with the results displayed in Figure 24. The highly tuned 
nature of any prototype design is exhibited here, where 
the input is of lower order than that for which the system 
was designed. An interesting result in Figure 24 which 
was not expected is the slightly decreased peak over- 
shoot in the case of the double-rate controller. 

The next design criterion applied was ripple-free 
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step response using both single and double-rate control- 
lers and the same error- sampling rate of the preceding 
design. The ripple-free process, K(z), was computed and 
supplied to the simulator, which calculated the single- 
rate controller, Dl. The double-rate controller was de- 
signed using the techniques discussed beginning at the 
bottom of page 20. Both controllers were quite simple 
and are given below: 

[)l - . Q 2 _ » Q QQ 4475 ~-- , 00000/ 43 d- 

I + 0-los- 3333 2”' 9 AO 

Figure 25 shows both systems responding without ripple to 
a noiseless step input. Again, the double-rate system is 
about one-half a second faster. The ramp response is next 
shown in Figure 26. Both systems have a steady state error 
since they were designed for a step input only. In the 
case of the double-rate controller the sustained ripple 
is due to the fact that the design technique was a special 
one, explicitly for ripple-free step response, and cannot 
be expected to give conventional behaviour for a different 
input. 

The results of the simulation studies using the second 
order approximation designs with inputs contaminated with 
Gaussian noise tend to reinforce the tentative conclusions 
mentioned earlier. It was impossible to observe the response 
of the continuous system to a noisy input since che contin- 
uous system was unstable. The disturbing result is that 
the sample-data systems, bcr,h with single and double-rate 
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digital controllers, which had given stable performance 
of varying degrees of goodness when the input was noise- 
free were, with but a single exception, unstable in the 
presence of noise. The exception is shown in Figure 27 
in which the design of D1 and D2 was for ripple-free re- 
sponse to a step and the input was a noisy step. It seems 
apparent at this point that a large signal-to-noise ratio 
is a sine qua non in the input to a sampled~data system. 

In the preceding experiments it was only the pulse 
transfer function of the plant and hold combination that 
was approximated in order to simplify the design calcula- 
tions. The differential equation of the plant itself 
remained of the fifth order for the numerical methods 
solution of all response values. 
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Figure 21 
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The simplicity of the design techniques for sampling rates 
of one or two samples per second which result in valid 
simplifying approximations is quite appealing, but there 
lurks an intuitive suspicion that these rates are too slow 
for a system required to accurately follow a high-speed 
aerial target. This suspicion has been confirmed in pri- 
vate communications between John 3. Slaughter and Gene F. 
Franklin, the results of which are that the damped natural 
frequency of the gun mount in train is about one and one-half 
cycles per second and that sampling should occur between six 
and twelve times the damped natural frequency of any system 
for good results. For these reasons a sampling rata of ten 
samples per second has been decided on by the group at N3L, 
and the remaining simulation studies of this work will use 
this rate. 

For T - 0.1 seconds the pulse transfer function of 
the plant and zero-order hold is shown below; 



Numerator of G(z) 


Denominator of 


21. 91 241 531 z" 1 


1.00000000 


131. 76730855Z -2 


-1. 61605197z -1 


63. S4130401z~ 3 


. 72871393z~^ 


3.03436539Z -4 


-.11377672Z -3 


.00419043Z -5 


,00612622z - ‘ i 




-. 00000152z"' J 



The roots of the numerator are located at the origin and at. 
-.00139903, -.05265621, -.47316625, and -5.48614151 

The presence of the zero of G(z) outside the unit 
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circle complicates the design effort slightly as it did 
in the third order plant considered earlier. The same 
methods apply and the resulting pulse transfer function 
of the process which gives prototype response to a step 
input is K(z) = .154174866Z" 1 + . 845825134z” 2 

Figure 28 shows this process responding to a noise- 
free step input, and Figure 29 shows the same process with 
a noise-free velocity input. The remarks made earlier con- 
cerning prototype processes apply and will not be repeated 
here. Figure 30 shows the results obtained with this de- 
sign when the step input contains additive noise. The 
best that can be said about the plant's behaviour in this 
case is that it remained stable and in the vicinity of the 
input signal. Not so for a noisy ramp input, to which the 
response of the sampled data system was unstable. It ap- 
pears that the sensitivity of a sampled data system to 
noise increases with the order of the input signal. 

The prototype design for zero steady state error in 
response to a velocity input begins with the equations below: 
K(z) = (1 + 5.4861415z" 1 )(a 1 z" 1 + a 2 z” 2 ) 

1 - K(z) = (1 - z _1 ) 2 (l + bjZ -1 ) 

From these are obtained three equations to be solved for 
three unknowns , namely 

a 1 = 2 - b 1# -5.4861415a 2 = b^ # and 
a 2 - 5. 486141 5a x = 2h ± -1 

Solving these, the resulting over-all process is 
K(z) = . 438754756z _1 + 2. 122490488z" 2 - 1. 561245244z~ 3 
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The coefficients of this process as well as the coeffi- 
cients of the numerator and denominator of G(z) and the 
value of the zero outside the unit circle were supplied 
to the simulator program, which then calculated a single- 
rate digital controller. The responses of the system to 
a noise-free ramp and a noise-free step are shown in Fig- 
ures 31 and 32 respectively. With the addition of noise, 
instability again resulted when the input was a rarnp, while 
for a noisy step the output was at least bounded for the 
duration of the experiment, although it began to look, in- 
auspicious towards the end. See Figure 33. In an attempt 
to achieve at least a semblance of stability in response 
to a noisy ramp signal the noise deviation was reduced 
from 0.2 to .05, but to no avail. A future investigation 
will be made into the stability threshold for noisy ramp 
signals. 

It was then decided to design for prototype step re- 
sponse using a staleness factor of 0.5. The design steps 
are shown below: 

Prototype K(z) = . 154174866Z - '*' + . 845825134z“ 2 as found 

earlier. This function contains the zero of G(z) located 

outside the unit circle at z = 5.48614151. The process 

pulse transfer function with a 0.5 staleness factor Is 

given by /-z >> _ ^2 (o^JSj4jn_4 & 2 ' + Q _ « 24_ 23_£_ i? 

~ / - a, sTH - ' 

For a unit step input, application of the Final Value Theo- 
rem shows that "a" must have the value 0.5. Thus K(z) is; 
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Kfc) = 



,^77^^33 Z -/ ^ 

/ - &. S’ Z~' 



Substituting in D/s) = QSjLl » — ^ (^ . ) . 

^ J P<z) I-KC2) 

where P(z) is the numerator of G(z) given on page 71, 



and 



Q(z) is the denominator of G(z) given on page 71, 

— 1 1 

and cancelling the common factors (1-z )(l+5.48614151z ), 

the single-rate digital controller is given by 



Numerator of D(z) 


Denominator of 


. 003517980 


1.00000000 


-. 002167290z _1 


. 95013414z -1 


, 000396307z -2 


. 24861928z -2 


-. 000021510z“ 3 


. 01088370z"” 3 


.000000025z"’ 4 


. 00001039Z -4 



This pulse transfer function was supplied to the simu- 
lator program? test inputs of noise-free step, noise-free 
ramp, noise-free acceleration, and noisy step were applied 
with the results shown in Figures 34 through 37 inclusive. 
The main point of interest in these results is that step 
response in the presence of noise for this staleness factor 
design is perhaps the best-behaved so far, although it is 
not good by any means. In addition, noisy inputs of ramp 
and acceleration were applied, but both of these gave un- 
stable responses as in previous designs. 

The final design effort to be described in this work 
was directed toward achieving ripple-free response to a 
velocity input when sampling at a rate of ten samples per 
second. For this design all of the zeros of G(z) must be 
contained in K(z), so we begin by writing 
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K(z) = ( l +5.48614151 z ~ 1 )( l +.47316625 z “ 1 ) (1+.05265&21 Z ” 1 ) 
(l+.ooisggosz " 1 )^^" 1 + a 2 z" 2 ) 

Since we are designing for zero steady state error for a 

ramp input, then it must be true that 

l-K(z) = (l-z" 1 ) 2 (l + b^ -1 + b 2 z“ 2 + b 3 z“ 3 + b 4 z“ 4 ) 

After multiplying out and equating like coefficients of z 
we obtain six simultaneous equation in six unknowns which 



are presented below in matrix form. 



1.0 


0 


1.0 


0 


0 


0 




a l 




2.0 


6.013363 


1.0 


-2.0 


1.0 


0 


0 




a 2 




- 1.0 


2.918063 


6.013363 


1.0 


-2.0 


1.0 


0 




b l 




0 


.140759 


2.918063 


0 


1.0 


-2.0 


1.0 




b 2 




0 


.000191 


.140759 


0 


0 


1.0 


-2.0 




b 3 




0 


0 


.000191 


0 


0 


0 


1.0 




b 4 




0 


— 
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This matrix equation was solved by computer using a Gaussian 
elimination library routine. With the results obtained there- 
from, the coefficients of the process, K(z), were calculated 

as follows: A(l) = 2 - = .319530814 

A( 2 ) = 2 b , - 1 - b 2 = 1.70120542 

A( 3 ) = 2b 2 - b 1 - b 3 = -.39202846 

A( 4 ) = 2b 3 - b 2 - b 4 = -.59772472 

A( 5 ) = 2b 4 - b 3 = -.03094094 

A( 6 ) = -b 4 = -.00004212 

These process coefficients were supplied to the simulator 
program along with the coefficients of P(z) and Q(z) and 
the values of the four plant zeros to be included. The 
program then calculated D(z) for the single-rate sampled 
data system. For a test input of a noise-free ramp the 
output has no overshoot or inter-sample ripple and reaches 
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zero error in 0.4 seconds. See Figure 38. The response 
to a noise-free step is shown in Figure 39 , where the 
"tuned" nature of the design is again exhibited. However, 
after the initial large overshoot the step response, too, 
is essentially ripple free and steady state in 0.4 seconds. 

For the noisy ramp signal shown in Figure 40 the sam- 
pled data system was unstable, while for a step signal with 
the same amount of additive noise the response can be seen 
in Figure 41. Note that the time scale in Figure 41 is five 
times longer than in the noise-free experiments, and the re- 
sponse to the noisy signal shows no signs of settling down 
any further at the end of five seconds than it did at the 
end of the first second. Clearly, noise of this power cannot 
be tolerated in the input signal to a sampled data system 
regardless of the digital controller design. 
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Figure 36 
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Figure 37 
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Figure 40 
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Figure 41 
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8. Conclusions. 

The behaviour of a unity feedback control system he. a 
been simulated using a digital computer program. U c mg this 
simulator program studies have been made of digital control 
methods for various plants, with particular attention fo- 
cused on the control of a Naval gun mount. The azimuth func- 
tion of the gun mount is described by a fifth order a infer- 
ential equation, and it is believed that the concic rons 
pertaining to it are applicable, at least in part, to other 
shipboard weapon systems. 

The use of multi-rate digital controllers wa- inves- 
tigated. In every'’ instance, the double-rate controller pro- 
vided faster response by one-half a sampling period. How- 
ever, only in the special case of ripple-free dec- . t.go for a 
step input was the double-rate response an unalloyed im- 
provement cvci the single-rate response. In genet 
transient concroj and inter-sample ripple was not uetice- 
ably improved using a double-rate controller fez the gun 
mount. For other plants of different frequency response, 
advantage mi gat oe taken of the fact that inter-sampie 
ripple when using a double-rate controller occurs ch t. ice 
the error- samp ling rate. 

Pulse transfer functions (or Z-transf orrns } ; a func- 

tion of the sampling period. It was found that i or certain 
sampling rates valid simplifying assumptions could made 
which facilitated the design procedure. Care should oe 
exercised that these simplifying assumptions are net made on 
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the basis of a sampling rate which is too slow in relation 
to the damped natural frequency of the system concerned. 

In this regard the approximations made to the pulse trans- 
fer function of the gun mount were inappropriate since they 
were based on a sampling rate roughly equal to the damped 
natural frequency instead of six to twelve times greater as 
must be the case for adequate control. 

Extreme sensitivity of sampled data systems to noise 
in the input signal was demonstrated. No satisfactory method 
of coping with input noise by design of the over-all process 
was found. So-called "variance reduction" processes showed 
little merit as the result of these simulation studies. The 
design criteria which came closest to producing satisfactory 
response to a noisy step input was the prototype with stale- 
ness factor of 0.5. 

Without compensation, either analog or digital, the 
gun mount is closed-loop unstable. Several digital control 
functions were developed which gave stable response to noise- 
free input signals. None of these gave stable operation when 
noise was present in a first order (ramp) input or a second 
order (acceleration) input, but for a noisy zero order (step) 
input digital control achieved a response that was at least 
bounded for the duration of the experiment, albeit unsatis- 
factory for most conceivable applications. It is concluded 
that sampled data system sensitivity to additive noise in the 
input increases with the order of the input signal. 

The computer program described in Appendix II for 
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finding Z-transf orms was very helpful in this investigation. 
It is hoped that this program, as well as the control system 
simulator upon the results of which this work was based, 
will prove useful to others in the field of sampled data 
systems. It is further hoped that this work has, in a small 
way, indicated the truthfulness of the assertion that digital 
compensation techniques introduce a degree of flexibility and 
accuracy in a control system that is not presently achievable 
with state-of-the-art analog techniques. 
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APPENDIX I 



DIGITAL COMPUTER PROGRAM FOR SIMULATION OF A SAMPLED 

DATA SYSTEM 



General Description of Program ALSTAPP 

This program is a simulation of the hybrid control 
system shown below. 




-3 



K(S) 






Three choices are available for the test input? namely , 
a unit step, a unit ramp, and a unit acceleration. The 
plant may be up to thirtieth order and is currently sub- 
ject to the following restrictions: 

(1) No plant zeros are allowed, i.e., the numerator 
of H(s) must be a constant. 

(2) The plant may not have a pole at the origin higher 
than first order, i.e., limited to Type 0 and Type 1 servos. 

(3) Only the constant term may be missing in the dif- 
ferential equation describing the plant. 

In essence, the plant is restricted to having a trans- 
fer function of the general form 






VK 



S(s+ ^5 + P 2 )*-« (5+ Pk-0 



simple poles. Pi 



H (s) = A(K+l) S K + A(k) s k_1 + • • • t >4(2) s +■ A (l) 

where A(l) is normally zero for Type 1 servo and is the only 

missing term allowed. 
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The rate at which the raw error, El, is sampled can 
be selected at will 0 This sample period T is chosen by 
the user and supplied to the computer. 

DELT is the small increment of time used in the simu- 
lation of the continuous plant by numerical solution of the 
plant's differential equation. This is done one thousand 
times, so that the total problem time is 1000 x DELT. 

DELT should be chosen as small as possible consistent with 
the desired total problem time. 

Although one thousand values of the solution are gen~ 
erated and are available for use by a computer plotting 
routine, only every tenth solution is produced on the Ana- 
lex printer, since it was felt that a column of one hundred 
values is sufficient for perusal by eye and more than ade- 
quate for plotting by hand. 

This program performs four major tasks. The first 
three are performed every time the program is executed, 
but it is optional with the user whether or not the fourth 
task will be done. These tasks are: 

(1) Simulation of a continuous, un-compensated, unity 
feedback system for the given plant. 

(2) Calculation of coefficients for a single-rate 
digital controller. 

(3) Simulation of the sampled-data system employing 
a single-rate digital controller using the coefficients 
calculated above. 

(4) Simulation, if desired, of a sampled-data system 
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employing a double-rate controller. The coefficients of 
the double-rate controller must be supplied by the user,, 
and if they are not, the fourth task is not performed. 

Information must be supplied to the program via data 
cards. Detailed instructions for the preparation and con- 
tents of these cards is given later. Seventeen (17) cards 
are required if the program is to do all four tasks* while 
only the first thirteen (13) need be supplied if simula- 
tion of the double rate controller is not desired. CAUTION § 
The cards must be presented to the computer in exactly the 
order and form specified. The slightest error will cause 
faulty operation or, perhaps, no operation at all. 

Discussion of Information Supplied by the User 

The general nature of the information that must be 
supplied to the computer is described next. 

(1) The user must choose the small time increment 
DELT to be used in the numerical method solution of the 
plant's differential equation. The smaller this increment 
the more accurate the solution. However, the total time 
simulated will be only one thousand times DELT, s© the user 
should first determine the time span of interest. . „ 1. e. , 

an estimate ©f how long it will take the system to achieve 
steady state, or how much of the transient condition it is 
desired to observe for lightly damped systems. . „ then divide 
this total time by 1000 to obtain the value for DELT. 

(2) The user must choose the sampling period to be 
employed in the sampled-data system with single-rate 
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controller simulation. This value is called T. 

(3) The user chooses the test input to be applied 

to the system. Unit values of step, ramp, and acceleration 
are available. 

{ NOTE Is These first three pieces of information, namely, 
DELT, T, and INPUT are supplied by data card number one. ) 

(4) The specifications of the continuous plant are 
supplied to the computer by data cards 2, 3, and 4. The 
order of the plant is denoted by K. This is the highest 
power of the Laplace variable, S, in the denominator of 
the plant transfer function, H(s). The numerator gain 
constant of H(s) is denoted by VK. The order of the pole 
at the origin is denoted by NTYPE and is limited to 0 or 1. 
These values, K, VK, and NTYPE are supplied by data card 
number two. 

(5) H(s) must be known in the polynomial form, that 

is to say, a numerator gain constant, VK, over a denominator 
polynomial in S, The constant coefficients of this denomi- 
nator polynomial are denoted by A(i) and must be supplied 
to the computer. The only coefficient which may have a 
zero value is the constant term A(l). The number of co- 
efficients is denoted by NA and is supplied by data card 
number 3. Even though A(l) is normally zero (in a Type 1 
servo) it must be included in the count, NA, on card 3 and 
its zero value shown in the first ten columns of data card 
number four (4). In the usual case of a iype 1 servo, the 
value of NA is one more than the order, K, of the plant. 
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(6) For the simulation of the sampled-data systems,, 

Z-transform techniques are implemented* Data reconstruc- 

tion is accomplished by a simulated zero order hold* which 

must be included in the Z-transform approach to the plant* 

If we call the plant and hold-circuit combination G(s), then 
- sT 

\ “ € and the pulse transfer func- 

s 

tion (or z-transform) is 



«... Z T ■ 8M 

This pulse transfer function for the plant and hold must 

be supplied by the user. Its general form is a ratio of 

-1 

polynomials in Z as shown below: 

“ MP 

G(z) = ■>■••• P(Nf) z 

’ 0(a) ~ Q(t) + +Q(3)£- 2 -+...^ 

where NP is the number of terms (and also the highest 
power of 7T ) in the numerator, P(Z), and NQ is the num- 
ber of terms in the denominator, Q(Z). In order for G(Z) 
to be a physically realizable pulse transfer function, the 
term Q(l) must be a constant, and in fact, is usually unity. 

(7) Although G(Z) is supplied to the computer in 
polynomial form, it must be known in factored form, in 
order that any poles or zeros of G(Z) that lie on or out- 
side the unit circle in the Z-plane may be accounted for. 
Usually, the plant is open loop stable and only zeros need 
be accounted for, but the program will also process infor- 
mation and provide the digital controller coefficients 
when the plant has a pole outside the unit circle. All 
such zeros must be supplied by the user under the generic 
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term ZERO(i) and the number of them is denoted by NZ. 

Up to ninety “-nine poles and/or zeros may be accounted for* 
(8) The overall pulse transfer function of the en- 
tire system is denoted by K(Z) and is defined as 

CO) 



K(z) * 



R(a) 



For the purposes of this program, K(Z) must be a 
numerator polynomial only. This implies a "finite-memory" 
process. If an infinite-memory or "non-finite settling 
time process" is desired, then it must be approximated 
by a finite number of terms of the infinite series ob- 
tained by long division. The program will operate on up 
to ninety-nine (99) such terms for K(Z). The general form 
for K(Z) is 

— I " N K 

«(2) - AK(i) 2 + AK(a) 2 ■+•***+ AK (nk) 2 

where NK is the number of terms provided. K(Z) must be 
supplied by the user. 

NOTE 2: In order to supply an array of numerical data to 

the computer two (2) data cards are required. The first 
card tells the computer how many values are in the array, 
and the next card (or cards) supplies these values. As- 
suming for the moment that no array contains more than 
eight (8) values, the following information is supplied 
to the computer by the data cards listed below s 

Card Nos. Information 

3 & 4 Array of coefficients of plant 

differential equation (i.e. , 
denominator of H( s ) ) . 



5 & 6 



Array of coefficients of nu- 
merator of G(Z) 0 
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Information 



Card Nos , 
7 & 8 



9 & 10 
11 & 12 



Array of values of zeros and 
poles of G(Z) to be included 
in overall process pulse trans- 
fer function, K(z), 

Array of coefficients of de= 
nominator of G(Z) e 

Array of coefficients of over= 
all process pulse transfer 
function, K(Z), 



(9) For certain investigations it may be desirable 
to limit the output of the digital controller t© a signal 
of a given magnitude to avoid plant saturation,, If so, 
the user must supply this value as a positive decimal 
value, E2KAX, and the program will limit the processed 
error from the controller to plus or minus this value. 

If no limiting is desired, a value of zero (0*0) should 
be given to E2I4AX and the program then sets a magnitude 
limit of one million on the processed error. Data Card 
13 contains the data on E2MAX. 

NOTE 3: In very simple cases it may be that there are no 

zeros or poles of G(Z) that need to be included in K(Z) 0 
When this occurs, obviously the datum NZ is zero and the 
array ZERO(i) does not exist. Under these conditions, 
prepare data card 7 as shown below 



C ? , l, 1 . 



Card 7 



00 



and remove card 8 completely! 
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Card No 



DATA CARDS for Input to Program ALSTAPP 



Format Contents Remarks 



1 


2F10 o 0 # Il 


DELT, T, INPUT 


Total time ■ 1000 
DELT 0 / Sample period 
T ss l y / error sampling 
rate© DELT is the in- 
crement used for simu- 
lation of continuous 
plant. INPUT * 1 for 
a step, 2 for a ramp, 

3 for acceleration 0 


2 


I2,F18.0,I1 


K, VK, NTYPE 


K = order of pi ant e 

VK - plant gain© 

NTYPE s type ©£ servo 

(io6 OJ 


3 


12 


NA 


NA “ no© terms in de- 
nominator of plant 
transfer function© 
Usually NA - K + 1. 

H(S) a 

VK 

A(K+1)S K + „. ©A(2)S+A(1) 


4 


8F10© 0 


A( I ) I = 1, NA 


A(I) is the array of 
coefficients ©f plant 
differential equation 
(i©e., denominator of 
H(S) above)© 


5 


12 


NP 


NP - no© terms in nu- 
merator ©f Z-transform 
of plant and zero- 
order hold whose z= 
transform is 

G(z) - PC.z Jt _ 

G{z} “ off! 


6 


8FX0©0 


P(I) 1=1, NP 


PU)z 4 +P(2)z" 2 t e . 

l+Q(2)z“ 1 "JQ(3)z“ 2 + BO © 

+ p(np)z~ np 

+Q(NQ) z-NQ+l 


7 


12 


NZ 


NZ = no© zeros of G(z) 
to be included in over- 
all pulse transfer fen, 
K(Z). 
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DATA CARDS for input to Program ALSTAPP ( Continued ) 



Card No. 


Format 


Contents 


Remarks 


8 


8F10.0 


ZERO(I) 1=1, NZ 


ZERO(I) = array of 
values of such zeros 
of G(Z) as are to be 
included in K(Z). 


9 


12 


NQ 


NQ = no. terms in de- 
nominator of G(Z). 


10 


8F10.0 


Q ( I ) 1=1, NQ 


Q( I ) = array of de- 
nominator coefficients. 


11 


12 


NK 


NK = no. terms in over- 
all pulse transfer 
function 

K(Z)=AK(l)z“ 1 +AK(2)z“ 2 
.. .+AK(NK)z“ NK 


12 


8F10.0 


AK( I ) I = 1, NK 


AK( I ) - array of pro- 
cess coefficients 
shown above. 


13 


F10.0 


E2MAX 


E2MAX is the limit 
placed on the absolute 
value of the output of 
the digital controller. 



Notes to the User : 

The FORMAT conveys information concerning the mode and 
location on the IBM card for each item of numerical infor- 
mation to be supplied to the program. If the following 
remarks are not sufficient, refer to the CDC Fortran Manual. 

"I" means an integer value. Whole numbers only. Do 
not show a decimal point. The integer value must be "right 
justified" in its allotted field. If the field is left 
blank, a value of zero will be inferred. 

"F" means a decimal value, either less than or greater 
than unity. A decimal point must be shown somewhere in the 
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allotted field*, If the value is negative, the minus sign 
must also fall within the field. Precise location of the 
value (including sign, if any, and the decimal point al - 
ways ) within the field does not matter. 

The numbers following I and F give the number of col- 
umns on an IBM card that are reserved for the field of 
that particular datum. For example, take the FORMAT con- 
trol for data card number 2: Columns 1 thru 2 are reserved 

for the integer value of K. If the integer value of K is 
four, then card 2 should have 04 in its first two columns. 
The next eighteen columns (3 thru 20) on card number 2 
must contain 1500000.0 somewhere between columns 3 and 20 
inclusive. On both cards 1 and 2, column 21 must contain 
a single integer. 

On those cards having an 12 format, columns 1 and 2 

contain an integer value between one and ninety-nine. If 
the value is less than ten, it must appear in column 2 
(i.e., the integer value must be right justified in its 
field. Decimal values need not be justified). In every 
case these 12 integer values tell the computer how many 
decimal values t© read from the next card or cards. This 
brings us t© the FORMAT 8F10.0. As before, the F1Q.0 
means a decimal value with ten card columns reserved for 
it. 8F10.0 means there may be up to eight such values on 
a single card (a card has eighty columns), and the com- 
puter will read another card if it has to, in order to 
read as many values as were specified by the previous 12 
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card. Most input arrays will contain eight or less values, 
however, and thus will fit nicely on a single card. Note 
that these decimal fields of ten columns each begin in 
card columns 1, 11, 21, 31, ...etc. Card 1 has two dec±~ 
mal fields of ten columns each, followed by an integer 
field of one (1) in column 21. Card 13 has a single 
decimal field occupying card columns 1 thru 10. 

If a decimal value of zero is to be entered, then 
place 0.0 in the appropriate field. For example, suppose 
array A has three decimal values A(l) - zero, A{2) * 
twenty, A(3) - minus two and a half. Then Card 4 should 
contain 

Between Columns The Decima l V alue 

1 and 10 0.0 

11 and 20 20.0 

21 and 30 -2.5 

If no sign is shown, the value is taken as being 

positive. 

In order to use this program to test a system using 
a single - rate controller, only the first thirteen data 
cards are required. 

If you desire to test a system using a doubl e- rate 
controller, then an additional four data cards are re- 
quired. These cards contain information as to the double- 
rate controller constants and follow card 13 (no inter- 
vening blanks 8 ) as shown below? 
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ADDITIONAL DATA CARDS FOR DOUBLE-RATE CONTROLLER 



Card No . Format Contents Sana rks 



14 


12 


NAA 


NAA is the no. terms 
in numerator of double- 
rate controller, D(Z^). 


15 


8F10.0 


AA( I ) 1=1 , NAA 


Array of numerator 
coefficients. 


16 


12 


NBB 


NBB is the no. terms 
in denominator of 

D(z 

2 BB(X)*BB{2)z“ 2 + 

.. .+aa(naa)z“ naa+1 
...SB(Na3)z- NBE+l 


17 


8F10.0 


BB( I ) 1=1, NBB 


Array of denominator 
coefficients. 



Example Problem usinq Program ALSTAPPg 

l | i item , ft — ■ I— ■ i fft I I -i ffi» i J ■ i ■i»A— — min im ■■ h i i r — i 



H(S) 



80, 00 0 

S 3 + 95S 2 + 1750S 



3y looking at H(s) we can see that VK = 80G00.0 

K ** 3 
NA * 4 
A(l) » 0.0 

A( 2 ) * 1750.0 
A( 3 ) * 95.0 
A(4) a 1.0 
NTYPE - 1 



Let us choose to observe the response of this plant for 
one second after the application of a unit step Input. 
So DELT ss 1/1000 - 0.001 
INPUT » 1 
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Let us select a sampling rate of 100 times per sec- 



ond for the sampled data system, i.e. , T * QoCl., 

We prepare the first four data cards as shown below* 



Column 

1 



Column 

11 



Column 

21 



31 




By partial fraction expansion and tables of Z-tr ans- 
forms we find the pulse transfer function of the plant- 

and-hold to ha 

in factored forms 

i i + 



G( Z ) = 

(1“S : 



209327Z ) (1 + 2 0 9 6727 3Z 

^0i778801Z~ 1 ) (1 



0 o 49658 5Z l ) 



in polynom Lai fcrm^ 

G(7) = « 010600Z°; 1 + „ 03 3672Z~ 2 + .0065 84Z" 3 . 

1.0 - 2.2753862“* + 1.662127Z* 2 - 0 o 3867412“ 3 

There is one aero outside the unit circle at -2*967275. 

Choosing a minimal prototype process we find that 

K(Z) “ 0*2520S2Z~ 1 + 0 o 747938Z“ 2 

and we see the following data: 

NP - 03 P(l) = .010600 

P( 2 ) = .033672 
P( 3 ) = .006584 

NZ s 0.1 ZERO (1) = -2.967275 
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NQ = 04 



NK - 02 
So we prepare 



Q(l) 
Q( 2) 
Q( 3) 
Q(4) 

AK(1) 

AK(2) 

data cards as 



= 1.0 

= -2.275386 
= 1.662127 
= -0.386741 

= 0.252062 
= 0.747938 

follows? 



Column Column Column Column 

1 11 21 31 




Card 12 
Card 11 
Card 10 
Card 9 
Card 8 
Card 7 
Card 6 
Card 5 



Suppose we do not desire to limit the output of the 

digital controller. Then set E2MAX = zero by preparing 

card 13 as shown below? 

Column 

1 

/ (0.0 Card 13 

( 

If it is not desired to simulate a double-rate digital 
controller * then only the first thirteen (13) data cards 
are required. 

If double-rate action is desired* then the user must 
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supply the coefficients of the double-rate controller. 

The pulse transfer function of the double-rate controller 
has the following general form: 

D(Z 0 ) = AA ( 1 * + AMSJZ" 1 + AA(3)Z“ 2 + ... + AA(NAA) Zj NAA+1 
BB( 1 ) + BB(2)Zj 1 + BBOJZ^ 2 + ... + B3(NBB)Z“ NBB+1 

where all coefficients have been normalized with respect 
to BB( 1 ) and BB(1) - 1.0. 

The format for cards 14 thru 17 is identical to that 
for cards 5 thru 12. 

In the double-rate controller the error-sampling 
period, T, remains the same as it was for the single-rate 
case. It is the output from the controller that occurs 
at the double-rate. 

Modification to allow for addition of white nois e to the 
test input . 

When used as described in the main body of this 
appendix, the program simulates a control system whose 
input is completely deterministic. Gaussian noise of 
selected mean and variance may be added to the input by 
making an addition to the first data card and inserting 
a new data card between the previous card 1 and 2 0 as 
follows. 

1. The presence of noise in the input must be re- 
quested by placing the integer 1 in column 31 of the first 
data card (which already contains DELT, T, and INPUT). 

2. The new data card, which must follow card 1 

above, contains an integer number between 01 and 99 right 
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justified in columns 1 and 2, and a decimal value in 
columns 3 through 20. The decimal value specifies the 
selected mean value of the noise, and in most cases will 
be shown as 0.0. 

The integer number specifies the number of uniform- 
ly distributed random numbers from which each normally 
distributed random number is obtained. For a good approx- 
imation to white noise, this number should be between ten 
and twenty. The variance of the Gaussian distribution 
depends on this number, N, as shown below? 

Variance - ^ 

As N approaches one, the noise approaches a uniform 
distribution between -1/2 and +1/2, for a selected mean of 

zero. 

The remainder of the data cards are then included 
in exactly the order and format described earlier. 

Modification to allow the program to use single - rate digi - 
tal controlle r coefficients supplied by the user . 

When used as described in the main body ©f this 
appendix, the program computes the coefficients of the 
single-rate controller using the following formula? 

0£Z) - . ! 2 kTz)' 

In order to do this, the user must supply data on 
the numerator and denominator of G(Z), the overall pulse 
transfer function, K{Z), and which zeros, if any, of G(Z) 
that are included in K(Z)„ CAUTION : Slide-rule accuracy 
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is usually not sufficient except in very simple cases i 
It is recommended that a desk calculator be used to pro- 
vide this data accurately to at least five decimal places. 

In case the user wishes to calculate his own con- 
troller coefficients for the single-rate sampled-data 
system simulation^ he may cause the program to omit the 
calculation of D(Z) and to employ the controller supplied 
by the user. In order to do this,, the data cards are ar- 
ranged as shown below? 



rd No. 


Contents 


Remarks 




1 


DELT, T, INPUT 


Same as before 




2 


K, VK 0 NTYPE 


Same as before 




3 


NA 


Same as before 




4 


A(I) I = 1, NA 


Same as before 




5 


Blank card! 


Indicates D(Z) supplied by 
user 


6 


NAA 


Single-rate controller 
efficients 


CO- 


7 


AA( I ) 1=1, NAA 


Single-rate controller 
efficients 


CO- 


8 


NBB 


Single-rate controller 
efficients 


co- 


9 


BB( I ) I - 1, NBB 


Single-rate controller 
efficients 


co- 


10 


32MAX 


Same as before 




11* 


NAA 


Double-rate controller 
efficients 


co- 


12* 


AA( I ) I ■- 1, NAA 


Double-rate controller 
efficients 


CO- 


13* 


NBB 


Double-rate controller 
efficients 


CO- 



Ill 



Card No. 



Contents 



Remarks 



14* 3B(I) I = 1, N3B Double-rate controller co- 

off icients 

♦Optional. 

The first four cards are just the same as they would 
be for a deterministic input where D(Z) is to be computed 
by the program. A blank fifth card signals the program 
that Z>(Z) will be supplied by the next four cards. These 
next four cards (6, 1, Q, and 9) give the data for tne 
single-rate controller exactly as described earlier in the 
main body of this memo for the double-rate controller co- 
efficients. Card 10 gives 3 2 MAX. If double-rate action 
is also desired, the last four cards must give the double- 
rate D(Z). 
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FLOW DIAGRAM for Program ALSTAP P 
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APPENDIX II 



DIGITAL COMPUTER PROGRAM FOR FINDING Z-TRANSFORMS AND 

ROOTS OF POLYNOMIALS 

General Description of Program Stappl 

This program is a combination rootfinder and Z-trans- 
form taker. It will operate in one of three modes , as 
specified by the user. These modes are: 

Mode 1 

Program operates as a rootfinder only. The user sup- 
plies the degree and coefficients of the polynomial whose 
roots are to be found. Handles polynomials up to and in- 
cluding thirtieth degree. 

Mode 2 

It is desired to take the Z-transform of a transfer 

function which is not known in factored form, such as 

H(s) = VK 

A(N+1)s N +A(N)s N " 1 + . . . + A( 2) s + A( 1 ) 

The program first finds the roots of the denominator 

polynomial. Then it takes the Z-transform of the above 

plant transfer function preceded by a zero-order hold, i.e. 

G < z > = & ( ¥' st - H(s)) - = §{§} 
and presents the z-plane poles of G(Z) as well as the 
numerator and denominator polynomials in Z”^. Then it 
operates as a rootfinder on the numerator, P(Z), in order 
to display the zeros of G(Z). 

Mode 3 

It is desired to find the Z-transform of a plant and 
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hold combination when the transfer function of the plant 
is known in factored form, such as 



As in Mode 2, the Z-transform is taken with the plant 
preceded by a zero-order hold, i.e. 



The program displays the Z-plane poles of G(Z) and the 
numerator and denominator polynomials. It then operates as 
a rootfinder on the numerator polynomial, P(Z), in order to 
display the zeros of G(Z). 

Restrictions . 

Mode 1 is relatively unrestricted. It will find the 
roots (real and imaginary parts) of polynomials in the 
following form (for N up to and including 30): 

P(x) = A(N+1) x N + A(H]x N= 1 + ... +A ( 2 ) x + A( 1 ) 

Modes 2 and 3 are restricted to polynomials of degree 
ten or less, with all roots lying on the real axis and no 
repeated roots. Mode 2 will accept a polynomial whose 
eleventh root is at the origin. In both modes the Z- 
transform operation is predicated upon the plant being a 
type 1 servo (i.e. the Laplace transfer function has a 
single pole at the origin). 

Since the program will not take the Z-transforrn 
when roots are complex, Mode 2 checks the magnitude of 
each imaginary part with respect to the corresponding real 
part of each root it finds. If the magnitude of any 
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imaginary part exceeds l/100th of the corresponding real 
part, the program stops. If all imaginary parts of the 
roots found in the first phase of Mode 2 are less than 
l/100th the magnitude of their corresponding real parts, 
then they are considered negligible and the root is treat- 
ed as entirely real. The program then proceeds as de- 
scribed earlier. 

If# thru error, Mode 2 is selected when the unfactored 
denominator of H(s) is higher than degree ten (exclusive 
of root at origin) then the program stops after finding 
the roots of the denominator. 

Mathematical Methods : Root - finder 

The rootfinding techniques is a modification of Bair- 

stow's iterative procedure for finding a quadratic divisor 
2 

of the form x -PP x-QQ. When this divisor is found suf- 
ficiently closely, it is solved by the quadratic formula, 
and divided out to leave a reduced polynomial two degrees 
lower. The criterion of suitable closeness is applied 
first to the change, DPP, in PP and then, if PP is suf- 
ficiently accurate, to the change, DQQ, in QQ D In both 
cases it is considered sufficiently accurate when the 
magnitude of the next correction to the coefficient is 
less than 10”^ of the magnitude of the coefficient itself. 

The essence of Bairstow's iterative procedure is to 
divide a polynomial, say A(x), by the trial quadratic 
divisor to obtain a first quotient polynomial, say Q(x), 
plus a first remainder. Then this first quotient is 
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divided again by the same trial divisor to obtain a sec- 
ond quotient, say T(x), plus a second remainder. Then 
the corrections necessary in the coefficients of the 
trial quadratic divisor are determined from a suitable 
combination of the coefficients of the first and second 
remainders. A symbolic example will illustrate the method. 

Let F(x) = A(N+1 )x N +A(N)x N “‘*‘+ . . . +A( 2 )x+A( 1 ) = the poly- 
nomial whose roots are to be found. 

2 

Let the trial quadratic divisor be TD(x) = x -PP.x-QQ 

NOW F(x) = TD(x)« rQ(N+lx N=1 +Q(N)x N " 2 + c . 0 +Q(4)x+Q(3)l 
+Q( 2 )x+Q( 1 ) L 

where the first quotient is Q(x) inside the brackets and 
the first remainder is Q(2)x+Q(l). 

Dividing again by the trial quadratic divisors 

Q(x) = TD(x)« [t(N-1)x N_3 +T(N-2)x N " 4 +...+T(4)x+T( 3)1 
+T( 2 )x+T( 1 ) 

where the second remainder is seen to be T(2)x+T(l) 

The corrections to the coefficients of TD(Z) are computed 
from the following formulae: 

M = PP.T(l) + QQ. T( 2 ) 

D = T(1 ) 2 -M. T( 2 ) 

Dpp = T( 2) «Q(1) - T( 1 ) *Q( 2) 

DQ q = M f Q(2) - T( 1 ) »Q.(1 | ) 

The modification of Bairstow°s procedure came about as 
the result of several runs using test polynomials. It 
was found that convergence in the case of non-repeated 
roots was more consistently obtained if each set of iter- 
ations on successive reduced polynomials was begun from 
the origin. In case there are repeated roots the return 
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to the origin is delayed until all roots at that location 
have been found. It was also found that an occasional 
set of iterations would not terminate even though no fur- 
ther refinements were being made in the coefficients of 
the trial divisor. When this occurs the program assumes 
the trial divisor is satisfactory after fifty iterations. 
No degradation of accuracy has been observed due to this 
artificial convergence device. 

Because of the sharpness of the criterion for a 
satisfactory trial divisor it appears that accuracy of 
root values to at least seven significant figures can be 
expected. Since the sum of the roots of an Nth degree 
polynomial is the coefficient of the term containing 
the variable to the (N-l)th power,, the program calcu- 
lates a check Siam from the root values found and dis- 
plays this check sum along with the true value as an in- 
dication of the accuracy achieved. This also provides 
the user with the means to tell at a glance whether or 
not the procedure converged. 

Mathematical Methods : Z-transformer 

The Z-transform portion of the program is restricted 
to type 1 servos preceded by a zero order hold. That is, 
the general foarm of the function whose Z-transform is to 
be taken must be 

-sT 

n , x _ 1-e S ' . VK 

“ S s( s-P x ) ( s-P 2 ) . . . ( s-P N ) 

Where s is the Laplace variable , T is the sampling period. 
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and N is the number of poles on the real axis? not in- 
cluding the pole at the origin. 

Note : 

All poles must be real and none can be repeated. 

sT 

By definition, Z=e , so by partial fraction expansion we 
can write 

G(Z) = VK>(Z-l)y/ A B Cyj. C£2l . C(N)\ 

Z A\s 2 S S - P l S = P 2 s “ Pn / 

The program first computes the partial fraction co- 
efficients, A, B, and C(i) i = 1, 2, . . . , N. The Z- 
transform of each term is then taken as follows: 




( eg) ) = c u) 

S - Pi Z - e PiT 

The poles of G(Z), exclusive of the one at unity, are 
PiT 

seen to be at Z = e for i = 1, 2, . .., N 

The program then proceeds to combine the individual 
Z-transforms over a common denominator, cancel the external 
factor (1-Z“^), and multiply through by the plant gain con- 
stant, VK. 

The final results are displayed as the coefficients 

of the numerator and denominator polynomials shown below: 

G( Z ) = P(PZ N + P(2)Z N ~ 1 + P(3)Z N ~ 2 + + P(N)Z + P(N+1) 

Q(1)Z N+1 + Q(2)Z N + Q( 3)Z N “ X + ... + Q(N+1)Z + Q(N+2) 

or, equivalently: 
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G(Z) = PU)^" 1 + P(2)Z~ 2 + P(N)Z“ N f P(N-H)z° ( ^ 1} 

Q(1)+Q(2)Z _1 +Q ( 3 ) Z” 2 + ... +Q(N+l)Z“ N +Q(N+2)Z” tN+1 ^ 



Discussion of the Data Provided by the User 

The general nature of the input data for each mode 
is discussed below. The sequence in which the data are 
presented is vital , and is correctly shown in the discussion. 
Mode Is 

1) MODE = the integer number 1. 

2) N = an integer value between 1 and. 30 representing 
the degree of the polynomial whose roots are to be found. 

3) A(i) = a linear array of decimal values repre- 
senting the coefficients of the polynomial whose roots are 
to be found. These values must appear on the data card(s) 
from left to right in the same order as they are in the 
polynomial when written in descending powers of the vari- 
able. Example for N = 4j 

A ( 5 ) x 4 + A(4)x 3 + A(3)x 2 + A(2)x + A(l) = F(x) 

Mode 2s 

1) MODE = the integer number 2. 

2) N = an integer value between 1 and 10 representing 
the degree of the denominator polynomial of the plant 
transfer function,, H(s). 

3) A(i) = a linear array of decimal values repre- 
senting the coefficients of the denominator of H(s) in 
descending powers of A. Example for N = 5s 

VK 



H(s) = 



A(6)s 5 + A(5)s 4 + A(4)s 3 + A( 3 ) s 2 + A(2)s + A(l) 
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For a type 1 servo, A(l) = 0.0 and N is allowed to 



be as high as 11. The pole at the origin will be removed 
by the program and N will be reduced by 1 when a zero 
value is read for A(l). It is optional with the user 
whether or not he removes the pole at the origin himself, 
or allows the program to do it. 

4) PERIOD and VK = both are decimal values. PERIOD 
is the reciprocal of the sampling rate for which the Z- 
transfer is to be taken. VK is the numerator gain con- 
stant of the plant transfer function, H(s). 

Mode 3j 

1) MODE = the integer number 3. 

2) N = the integer number of poles of the plant 
transfer function, exclusive of the pole at the origin. 

3) P(i) = a linear array of the decimal values of 

the poles of H(s). Does not include the pole at the origin. 
All poles must be simple and must lie on the real axis of 
the s-plane. 

4) PERIOD and VK = same as in Mode 2 discussed above. 

Detailed Instructions for Preparing Input Data 

All input data are provided to the program through 
the use of IBM cards. The "normal" number of data cards 
required by each mode of operation is tabulated below: 



Mode 

1 

2 

3 



No . of data cards required 

3 

4 
4 
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r _ 1 

i 









"Normally" it is expected that the number of decimal 
values in the linear array will be eight (8) or less, so 
that the entire array can be contained on a single card 0 
For arrays of more than eight values, additional cards 
are used as necessary. 

For convenience, the following table summarizes the 
requirements for input data: 



TABLE A 



Card No. 


Datum 


Format 


Used in 


Mode 1 


Mode 2 


Mode 3 


1 


MODE 


11 


X 


X 


X 


2 


N 


12 


X 


X 


X 


3* 


( Array ) 


8F10.0 


X 


X 


X 


4 


PERIOD, VK 


F10.0, F20.0 




X 


X 



* Use additional cards 3 as necessary for arrays contain- 
ing more than eight (8) values. 

The "Format" information needs to be explained. The 
format letter I means that the value is an integer. In- 
tegers must be shown without any decimal point and must 
be "right justified" in their allotted space. Right 
justified means that the right-most digit of the integer 
value must appear in the right-most column of the space 
allotted for that value on the card. This applies only 
to integer values. 

The format letter F indicates a decimal value and a 
decimal point must be used . If the value is negative, 
the minus sign must appear. Unsigned values are taken 
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to be positive. The decimal value,, including the deci- 
mal point always and minus sign if necessary, may appear 
anywhere within the allotted space (or "field 1 "). 

The number following the letter I or F is the number 
of columns on the IBM card forming the allotted space or 
field for that value. There are eighty (80) columns on 
a card and they are assigned from left to right in order. 

When a number precedes the format letter, it tells 
how many times that particular field is repeated. Thus 
card 3 contains eight identical fields of ten columns 
each. The next card may be read by the program as a con- 
tinuation of card 3 if more than eight values are in the 
array. 

Card 1 uses only the first column. This column must 
contain either the digit 1, 2, or 3. 

Card 2 uses only the first and second columns. If 
the value of N is nine or less, then skip column one and 
place the digit in column two. 

Card 3 uses as many identical fields of ten columns 
each as are necessary. Note that these fields begin in 
columns 1, 11, 21, 31, etc. 

Card 4 uses the first ten columns to contain the deci- 
mal value for the PERIOD. The next twenty columns (i.e. , 

11 through 30) are reserved for the decimal value of VK. 
Decimal values may lie anywhere within their allotted 
field as long as the decimal point is shown in the field 
and the minus sign is in the field for negative values. 
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Let us now show the preparation of data cards for 
some representative problems. 

EXAMPLE PROBLEM - Mode 1 

Find all the roots of 

F(x) = x 7 + 16x 6 + 102x 5 + 371x 4 + 838x 3 + 1200x 2 + 1012x+ 420 
Since all we need is a rootfinding operation, we choose 
Mode 1. 

The degree, N, of the polynomial is 7. The polyno- 
mial is written in descending powers of x so the three 
data cards would be 



Column 


Column 


Column 


Column 


Column 


Column 


Column 


Column 




1 


11 


21 


31 


41 


51 


61 


71 




'1.0 

i 


16.0 


102.0 


371.0 


838.0 


1200.0 


1012.0 


420.0 ( 


Card 
1 3 


J 07 

I 
















, Card 
) 2 


H 


, 














[Card 



1 



The answers are: 

Roots 1 and 2 = -1.0 + j 1.0 

3 and 4 = -1.5 + j 1.658312 

5 and 6 = -2.0 + j 1.414214 

7 = -7.0 



EXAMPLE PROBLEM - Mode 2 



Find the Z-transform of the following configuration 

zero-order 

hold plant 



O 

20 

where the 



>1 




samples/sec 20 samples/sec 

plant transfer function is 
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H( s ) = — = 1,6064235 x_10^ 

s + 134s + 5064s J + 73274s^ + 356983s 

Since H(s) is not known in factored form we choose 

Mode 2. 

The degree, N, of the denominator of H(s) is 5*. 

The denominator of H(s) is written in descending 
powers. Note that A(l) is missing! 

The numerator gain constant, VK, is 1.6064235 x 10 9 . 
The sampling rate is twenty samples per second so 
the sampling PERIOD is .05 seconds. 

The data cards are prepared as shown below? 



Column 


Column 


Column 


Column 


Column Column 


Column 


Column 




1 


11 


21 


31 


41 51 


61 


71 




* 

1 .05 


1606423500.0 










Card 

4 


*1.0 

i 


134.0 


5064.0 


73274. 


0 356983.0 0.0 






Card 

3 


1 05 
1 














Card 

2 


1 2 
i 














Card 

1 



f 



* User may remove pole at origin, if desired, and enter a 

4th degree denominator where A(l) = 356983. 

The answers are? 

z-transform = o(z) = ^ 

Q(Z) Q(1)+Q(2)Z" 1 +Q(3)Z“^+Q(4)Z“ J + 

P(5)Z~ 5 

Q(5)Z” 4 + Q(6)Z" 5 
where 

P(l) = 1.54749758 Q(l) = 1.0 

P( 2 ) = 16.74877866 Q(2) = -2.33676594 
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P(3) = 16.29465286 Q(3) = 1.92220653 

P( 4 ) = 2.04809669 Q(4) = - .67242310 

P(5) = .01866472 Q(5) = .08821342 

Q( 6 ) = - .00123091 

The zeros of G(Z) are at Z = - .00988230 

- .13606251 

- .91925736 
-9.75793506 

The poles of G(Z) are at Z = 0.57694981 

0.42741498 

0.31663673 

0.01576442 

EXAMPLE PR03LSM - Mode 3 

Find the Z-transform of the following configurations 



zero -order 

hold plant 




where the plant transfer function is 

_ 1.6064235 x 10 9 

' ' s(s + 11) (s + 17) (s + 23) (s + 83) 

Since H(s) is known in factored form we choose Mode 3. 

The number, N, of real poles (not counting the one at 

the origin) is 4. 

g 

The numerator gain constant, VK, is 1.6064235 x 10 . 

The sample PERIOD is .05 seconds. 

The data cards are prepared as shown below: 

Column Column Column Column Column Column Column Column 

1 11 21 31 41 51 61 71 



i 

1 .05 

I 



1606423500.0 
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Card 

4 



Column Column Column Column Column Column Column Column 



1 11 21 31 41 51 61 71 

I — 

1-11.0 -17.0 -23.0 -83.0 Card 

! 3 

1 04 Card 

' 2 

I 1 1 

, 3 Card 

i 1 



The answers are the same as those obtained in the example 
problem for Mode 2. 
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FLOW DIAGRAM FOR PROGRAM STAPP1 
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APPENDIX III 



ftg'tlon of stability of a third -order 
System as s ampling rate is decreased 






or 

P(z) 

Q(») 






cr 

s + y+jj- + ” 5 



2-1 

+ 70 J 



Ala 4. Bz 

Fo* — 



a - 1 



Ca , 

- jrr + 

a-e 2 -e 



where K- g'O, 000 

As SI. / 4 j?X/ 0 ~ S ~ 
B- - 3. /OJO * /o~ r 

2 s_T c = ^.s-ss-x/tr* 

J £> = - 0 . 4 O 3 S*l 0 ~^ 



Mult out : 

A/ . 3 r -jrr -7or _ i r -arr -7cr -nri 

Q( 2 ) - 2 - |_<£ + G flj 2 + [s +e ■+■ <£ J 2 " — s 



-?rr 



-7oTi 



?(b) = k[at- Q(i as le c l i)- c( 2 + e OT J - P(^+e arT J] ^ 

+ K [- AT (e*l e ° T ) t B(e 3f +- e 7o7 i- i’ rr ) + C ( I + ;?? dT ) +• D(i+ 3 e**)] Z 

+ K[ATe 9rr - Be W - Ce 7 ‘ T - Pi^ T ] 






T*(t) + 



T 



o- 



l-£ rs 




£ 0 , 000 


s 




5 ’( 5 >;?.rJ('s> 7 o y ) 



c(-fc) 

— 



Characteristic equation: A( 2 )- ft*) + Q( 2 ) = £3 2 + G.& + 

where <£? 3 = 1 

^ = A'Ar- K" <3 - KB ~ OT - JKC-KC Q°- 3 KD - KD~& T - KB - e^-^° T - 1 

A 1 = -KATe*- KATe° T + KBz sr + KB£ 7or + Kb£ f ’ rT i- KC + JKc£~ toT 
+ kd -r iKDe s * r +. e jfr + « 7or + e ^ r 

^o = KAT e 9r - KB - KC e 7 °- KD^ rr - T 
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Note that G(z) has been expressed in positive powers of 
z since we anticipate making a change of variables. Later 
on„ z~transforms are expressed using inverse powers of z 
in order to correspond to the concept of difference equa= 
tions. 

_ £ -h u> 

Using the Bilinear Transformation? z — j~Z~Zo 

aw - &»(■££? + ^ 

— $ j ■*- to J «/- [u 40 + °° J f / ~ 40 ~ ^ + 10 ] p ~ ~~ 

= f to 2 [34 5 -4t-4+34*]+<o[3a s + 34 0 ] + 

Routh-Hurwitz criterion for third order equations 

all roots in left half plane if \^0 2 — b c b 3 > O 
to^re - &?> 4- +* &1 +• a o 

*b, = 3^<3 4* 0 .% & i “* 3^0 

b a = 3^3 -tfa + 3£* 

b 3 = Cl’l * H- j — ^ o 
X'K ter>vt 5 <3-f tlve 5yjte>^ |5Ara-weter5 3 

fe.. € W V^r; + i 7 ° T (-KAT) ♦ £ ST (-MT) i- (MT-KC-KO) 

b t = & r ‘(-3KfiT + 2K& + 2.) + e (V^T - J KO -3 ) 

t <£ arT (l<AT- 5KS- 3) 4- (KAT -3KC - 3KD -t a) 

s ^*'(3 KAT - 4 K&- A) + € (kAT- 4 *< c ) 

4- e 2rT ( KAT - 4 K D) 4- (-KAT tKC+KD + l) 
b 3 = e 9 ^ T (KAT) ^'g T (-KAT UK8t2KC4-2) 

4- e aST (-KAT 4 3K6 + 2KO 4- 2.) + (- KAT +-3KC 4 3KD 4 1 ) 

A computer program to solve the above equations for 
one hundred values of T and the results obtained are shown 
on the following pages. 
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. .JOB STAPP2 

PROGRAM STAPP2 

C DETERMINATION OF ROUTH-HURWITZ CRITERION FOR THIRD 

ORDER SYSTEM 

C AS SAMPLING PERIOD IS VARIED 

VK=80000. 

A= . 000571428 
B=-. 000031020 
C- .000035555 
D=-. 000004535 
T=0 . 0 
N~0 

PRINT 199 

199 FORMAT( 32H PERIOD ROUTH CRITERION) 

100 T~T+. 001 

EXP25s2.71828**(~25.0*T) 

EXP70— 2 , 71828** (-70. 0*T) 

EXP95=2.71828**(-95.0*T) 

AO=EXP95*(VK*A*T-VK*B-l. )-EXP70*VK*C-EXP25*VK*D 
0Al=EXP95*( VK*B+1. )+EXP70*(~VK*A*T+VK*B+2. *VK*C+1. ) 
1 EXP25*(~VK*A*T+VK*B+2.0*VK*D+1. )+VK*C+VK+D 

0A2=EXP70*( -VK*B-VK*C“1. )+EXP25*(-VK*B~VK*D~l. )+VK* 
(A*T“2. *C~2. *D) 

1 -1.0-VK*B 

A3— 1. 0 

B0=A3+A2+Al=A0 

Bl= 3 . 0 *A3+A2-Al- 3 . 0 * AO 

B2=3„0*A3‘~A2~Al+3.0*A0 

B3-A3-A2+Al“A0 

ROUTH* Bl*B2~B0*B3 

PRINT 200, T, ROUTH 

200 FORMAT( F10. 4, F20.5) 

N=N+1 

IF(N~100 )100, 300, 300 
300 STOP 
END 
END 



(FORTRAN program for the CDC 1604) 



PERIOD 


ROUTH CRITERION 


PERIOD 


ROUTH CRITI 


.0010 


.00060 


.0510 


“2. 33345 


.0020 


.00418 


.0520 


=2.44305 


.0030 


.01226 


.0530 


-2.55283 


.0040 


.02524 


.0540 


“2.66273 


.0050 


.04279 


.0550 


“2.77266 


.0060 


.06411 


.0560 


“2.88256 


.0070 


.08814 


.0570 


“2.99235 


. 0080 


.11372 


.0580 


“3.10199 


.0090 


.13966 


.0590 


=3.21141 


.0100 


.16483 


.0600 


“3.32057 


.0110 


.18818 


.0610 


“3.42940 


.0120 


.20879 


.0620 


“3. 53787 


.0130 


.22587 


.0630 


“3.64593 


.0140 


.23874 


.0640 


=3.75354 


.0150 


.24688 


.0650 


“3.86066 


.0160 


.24986 


.0660 


“3.96725 


.0170 


.24736 


.0670 


“4.07329 


.0180 


.23919 


.0680 


“4.17874 


.0190 


.22522 


.0690 


“4.28358 


.0200 


.20541 


.0700 


“4.38776 


.0210 


.17977 


.0710 


“4.49128 


.0220 


.14839 


.0720 


=4.59410 


.0230 


.11137 


.0730 


“4.69620 


.0240 


.06889 


.0740 


“4.79756 


— ► .0250 


STABILITY. 02112 


.0750 


“4.89816 


.0260 


LIMIT ~. 03171 


.0760 


=4.99798 


.0270 


-.08939 


.0770 


-5.09700 


.0280 


-.15167 


.0780 


“5.19522 


.0290 


-.21831 


.0790 


“5.29261 


.0300 


=. 28905 


.0800 


“5. 38915 


.0310 


-.36365 


.0810 


“5.48485 


.0320 


-.44185 


.0820 


“5.57967 


.0330 


=.52341 


.0830 


“5.67363 


.0340 


-.60809 


.0840 


-5.76669 


.0350 


-.69565 


.0850 


“5.85886 


.0360 


“.78588 


.0860 


“5.95012 


.0370 


-.87855 


.0870 


=6.04048 


.0380 


-.97345 


.0880 


“6.12991 


.0390 


“1.07039 


.0890 


“6.21842 


.0400 


“1.16919 


.0900 


=6.30599 


.0410 


-1.26965 


.0910 


“6.39264 


.0420 


“1.37162 


.0920 


“6.47834 


.0430 


“1.47493 


.0930 


“6.56309 


.0440 


“1.57944 


.0940 


“6.64690 


.0450 


-1.68499 


.0950 


“6.72977 


.0460 


=1.79146 


.0960 


“6.81168 


.0470 


“1.89873 


.0970 


“6.89263 


.0480 


“2.00666 


.0980 


“6.97264 


.0490 


=2.11516 


.0990 


“7.05169 


.0500 


-2. 22412 


.1000 


“7.12979 



Determination of Routh-Hurwitz criterion for a third 
order system as sampling period* T* goes from .001 to 0.1 
seconds. 
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